REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  NO.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments 

regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggesstions  for  reducing  this  burden,  to  Washington 

Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington  VA,  22202-4302. 
Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  oenalty  for  failing  to  comply  with  a  collection  o1 
information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


2.  REPORT  TYPE 
Final  Report 


1.  REPORT  DATE  (DD-MM-YYYY) 
10-02-2012 


4.  TITLE  AND  SUBTITLE 

COMPUTATIONALLY  EFFICIENT  MODELING  OF 
HYDROCARBON  OXIDATION  CHEMISTRY  AND  FLAMES 
USING  CONSTITUENTS  AND  SPECIES 


6.  AUTHORS 
Josette  Bellan 


3.  DATES  COVERED  (From  -  To) 

19-May-2008  -  30-Nov-2011 


5a.  CONTRACT  NUMBER 
W911NF-08-1-0130 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 
611102 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAMES  AND  ADDRESSES 

California  Institute  of  Technology 

Sponsored  Research  MC  201-15 

California  Institute  of  Technology 

Pasadena,  C A  91125  - 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND 
ADDRESS(ES) 

U.S.  Army  Research  Office 
P.O.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 
ARO 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

54125-EG.ll 


12.  DISTRIBUTION  AVAILIBILITY  STATEMENT 
Approved  for  Public  Release;  Distribution  Unlimited 


13.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  contrued  as  an  official  Department 
of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


14.  ABSTRACT 

This  report  describes  a  study  performed  under  ARO  sponsorship,  addressing  the  investigation  of  a  novel  way  to 
reduce  complex  and  extensive  oxidation  reaction  mechanisms  for  fuel  mixtures  containing  hundreds  of  species  to  a 
much  smaller  number  of  progress  variables,  typically  by  a  factor  of  ten.  The  study  has  also  been  extended  to 
computing  laminar  flames.  Because  the  results  have  been  documented  in  several  papers  published  in  the  refereed 
literature  and  manuscripts,  this  final  report  is  in  the  form  of  an  Executive  Summary  succinctly  describing  the  results 


15.  SUBJECT  TERMS 

chemical  kinetics  mechanism  reduction;  laminar  flame  prediction 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

15.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

h.  ABSTRACT 

c.  THIS  PAGE 

ABSTRACT 

OF  PAGES 

Josette  Bellan 

UU 

UU 

UU 

UU 

19h.  TELEPHONE  NUMBER 

818-354-6959 

Standard  Form  298  (Rev  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 


Report  Title 

COMPUTATIONALLY  EFFICIENT  MODELING  OF  HYDROCARBON  OXIDATION  CHEMISTRY  AND 
FLAMES  USING  CONSTITUENTS  AND  SPECIES 


ABSTRACT 

This  report  describes  a  study  performed  under  ARO  sponsorship,  addressing  the  investigation  of  a  novel  way  to  reduce  complex  and 
extensive  oxidation  reaction  mechanisms  for  fuel  mixtures  containing  hundreds  of  species  to  a  much  smaller  number  of  progress  variables, 
typically  by  a  factor  of  ten.  The  study  has  also  been  extended  to  computing  laminar  flames.  Because  the  results  have  been  documented  in 
several  papers  published  in  the  refereed  literature  and  manuscripts,  this  final  report  is  in  the  form  of  an  Executive  Summary  succinctly 
describing  the  results  and  putting  them  in  perspective  with  respect  to  existing  oxidation  mechanism  reduction  schemes  and  flame 
computations.  The  papers  and  manuscripts  are  individually  listed  as  Appendices,  and  attached  to  this  report. 


Enter  List  of  papers  submitted  or  published  that  aeknowledge  ARO  support  from  the  start  of 
the  projeet  to  the  date  of  this  printing.  List  the  papers,  ineluding  journal  referenees,  in  the 
following  eategories: 

(a)  Papers  published  in  peer-reviewed  journals  (N/A  for  none) 


Received 
2012/02/10  II  10 


Paper 

Kenneth  Harstad,  Josette  Bellan.  A  model  of  reduced  kinetics  for  alkane  oxidation  using  constituents  and 
species:  Proof  of  concept  for  n-heptane,  Combustion  and  Flame,  (08  2010):  0.  doi: 

1 0.1 01 6/j.combustflame. 201 0.02.01 3 


2010/10/04  01  3  K.  Harstad,  J.  Bellan.  A  model  of  reduced  oxidation  kinetics  using  constituents  and  species:  iso-octane  and  its 
mixtures  with  n-pentane,  iso-hexane  and  n-heptane.  Combustion  and  Flame,  (03  2010):  .  doi: 


TOTAL:  2 

Number  of  Papers  published  in  peer-reviewed  journals: 


(b)  Papers  published  in  non-peer-reviewed  journals  (N/A  for  none) 


Received  Paper 


TOTAL: 

Number  of  Papers  published  in  non  peer-reviewed  journals: 


(e)  Presentations 

“Alkane  Kinetics  Reduction  Consistent  with  Turbulence  Modeling  using  Large  Eddy  Simulation”,  (K.  G.  Harstad  and  J.  Bellan), 
AIAA-2010-1514,  presented  at  the  48th  Aerospace  Sciences  Meeting,  Orlando,  FL,  January  4-7,  2010 

“Computation  of  Laminar  Premixed  Flames  Using  Reduced  Kinetics  Based  on  Constituents  and  Species”,  (K.  G.  Harstad  and  J.  Bellan), 
AIAA-201 1-ID892667,  presented  at  the  49th  Aerospace  Sciences  Meeting,  Orlando,  FL,  January  4-7,  2011;  also  paper  1E08  presented  at 
the  7th  US  National  Combustion  meeting,  Atlanta,  GA.,  March  21-23,  201 1 

“Modeling  of  Steady  Laminar  Flames  for  One-dimensional  Premixed  Jets  of  Heptane/Air  and  Octane/Air  Mixtures”,  (K.  G.  Harstad  and  J. 
Bellan),  AIAA-20 12-0340,  presented  at  the  50th  Aerospace  Sciences  Meeting,  Nashville,  TN,  January  9-12,  2012 

Number  of  Presentations :  3.00 


Non  Peer-Reviewed  Conferenee  Proeeeding  publieations  (other  than  abstraets): 


Received  Paper 

2012/02/10  1'  9  K.  G.  Harstad,  J.  Bellan.  Modeling  of  Steady  Laminar  Flames  for  One-dimensional  Premixed  Jets  of 
Heptane/Air  and  Octane/Air  Mixtures,  Aerospace  Sciences  Meeting.  2012/01/09  03:00:00,  .  :  , 

2012/02/10  i:  8  J.  Bellan,  K.  G.  Harstad.  Alkane  Kinetics  Reduction  Consistent  with  Turbulence  Modeling  using  Large  Eddy 
Simulation,  Aerospace  Sciences  Meeting.  2010/01/09  03:00:00,  .  : , 

201 1/08/22  11  6  Kenneth  G.  Harstad,  Josette  Bellan.  Computation  of  Laminar  Premixed  Flames  Using  Reduced  Kinetics  Based 
on  Constituents  and  Species,  Aerospace  Sciences  Meeting.  2011/01/04  03:00:00,  . :  , 

TOTAL:  3 

Number  of  Non  Peer-Reviewed  Conference  Proceeding  pubiications  (other  than  abstracts): 

Peer-Reviewed  Conferenee  Proeeeding  publieations  (other  than  abstraets): 


Received  Paper 

TOTAL: 

Number  of  Peer-Reviewed  Conference  Proceeding  pubiications  (other  than  abstracts): 


(d)  Manuscripts 

Received  Paper 

2011/08/22  1‘  5  Kenneth  G.  Harstad,  Josette  Bellan.  A  model  of  reduced  oxidation  kinetics  using  constituents  and  species: 

iso-octaneand  its  mixtures  with  n-pentane,  iso-hexane  and  n-heptane.  Combustion  and  Flame  (03  2010) 

2010/12/20  01  4  K.  Harstad,  J.  Bellan.  Computation  of  Laminar  Premixed  Flames  Using  Reduced  Kinetics  Based  on 
Constituents  and  Species,  (12  2010) 

2010/03/31  2  K.  Harstad,  J.  Bellan.  A  model  of  reduced  oxidation  kinetics  using  constituents  and  species:  iso-octane  and  its 

mixtures  with  n-pentane,  iso-hexane  and  n-heptane,  (03  2010) 

2010/03/31  1‘  1  K.  Harstad,  J.  Bellan.  A  model  of  reduced  kinetics  for  alkane  oxidation  using  constituents  and  species:  proof  of 
concept  for  n-heptane  ,  (03  201 0) 

TOTAL:  4 

Number  of  Manuscripts: 


Books 


Received 

TOTAL: 


Patents  Submitted 


Patents  Awarded 


Awards 


None. 


Graduate  Students 


NAME  PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 


Names  of  Post  Doctorates 

NAME 

PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 

Names  of  Faculty  Supported 


NAME 

FTE  Equivalent: 

Total  Number: 

PERCENT  SUPPORTED 

Names  of  Under  Gradnate  stndents  supported 

NAME 

PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 

Student  Metrics 

This  section  only  applies  to  graduating  undergraduates  supported  by  this  agreement  in  this  reporting  period 

The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period: .  0.00 

The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period  with  a  degree  in 

science,  mathematics,  engineering,  or  technology  fields: .  0  00 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  continue 

to  pursue  a  graduate  or  Ph.D.  degree  in  science,  mathematics,  engineering,  or  technology  fields: .  0.00 

Number  of  graduating  undergraduates  who  achieved  a  3.5  GPA  to  4.0  (4.0  max  scale): .  0.00 

Number  of  graduating  undergraduates  funded  by  a  DoD  funded  Center  of  Excellence  grant  for 

Education,  Research  and  Engineering: .  o.OO 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  intend  to 

work  for  the  Department  of  Defense .  0.00 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  receive 

scholarships  or  fellowships  for  further  studies  in  science,  mathematics,  engineering  or  technology  fields: .  0.00 

Names  of  Personnel  receiving  masters  degrees 

NAME 

Total  Number: 


Names  of  personnel  reeeiving  PHDs 

NAME 

Total  Number: 


Names  of  other  researeh  staff 


NAME 

PERCENT  SUPPORTED 

Josette  Belan 

0.05 

K.  G.  Harstad 

0.35 

FTE  Equivalent: 

0.40 

Total  Number: 

2 

Sub  Contractors  (DD882) 


Inventions  (DD882) 


Seientifie  Progress 


Described  in  the  attachment. 


Technology  Transfer 


FINAL  REPORT 


COMPUTATIONALLY  EFFICIENT  MODELING  OF  HYDROCARBON 
OXIDATION  CHEMISTRY  AND  FLAMES  USING  CONSTITUENTS  AND 

SPECIES 


Josette  Bellan 

Mechanical  Engineering  Department 
California  Institute  of  Technology 
Pasadena  CA  91125 


Abstract 


This  report  describes  a  study  performed  under  ARO  sponsorship,  addressing  the  investigation 
of  a  novel  way  to  reduce  complex  and  extensive  oxidation  reaction  mechanisms  for  fuel  mixtures 
containing  hundreds  of  species  to  a  much  smaller  number  of  progress  variables,  typically  by  a 
factor  of  ten.  The  study  has  also  been  extended  to  computing  laminar  flames.  Because  the  results 
have  been  documented  in  several  papers  published  in  the  refereed  literature  and  manuscripts,  this 
final  report  is  in  the  form  of  an  Executive  Summary  succinctly  describing  the  results  and  putting 
them  in  perspective  with  respect  to  existing  oxidation  mechanism  reduction  schemes  and  flame 
computations.  The  papers  and  manuscripts  are  individually  listed  as  Appendices,  and  attached  to 
this  report. 


TABLE  OF  CONTENTS 


EXECUTIVE  SUMMARY . 1 

REFERENCES . 4 

APPENDICES . 5 

Appendix  1:  . 5 


“A  model  of  reduced  kinetics  for  alkane  oxidation  using  constituents  and  species:  proof 
of  concept  for  n-heptane”,  (K.  G.  Harstad  and  J.  Bellan),  Combustion  and  Flame,  157, 
1594-1609,  2010 

Appendix  2: . 6 

“A  model  of  reduced  oxidation  kinetics  using  constituents  and  species:  iso-octane 
and  its  mixtures  with  n-pentane,  iso-hexane  and  n-heptane”,  (K.  G.  Harstad  and  J. 
Bellan),  Combustion  and  Flame,  157,  2184-2197,  2010 

Appendix  3: . 7 

“Alkane  Kinetics  Reduction  Consistent  with  Turbulence  Modeling  using  Large  Eddy 
Simulation”,  (K.  G.  Harstad  and  J.  Bellan),  AIAA-2010-1514,  presented  at  the  48**^ 
Aerospace  Sciences  Meeting,  Orlando,  FL,  January  4-7,  2010 

Appendix  4: . 8 

“Computation  of  Laminar  Premixed  Flames  Using  Reduced  Kinetics  Based  on 
Constituents  and  Species”,  (K.  G.  Harstad  and  J.  Bellan),  AIAA-201 1-ID892667, 
presented  at  the  49*  Aerospace  Sciences  Meeting,  Orlando,  FL,  January  4-7,  2011;  also 
paper  1E08  presented  at  the  7*  US  National  Combustion  meeting,  Atlanta,  GA.,  March 
21-23,2011 

Appendix  5: . 9 

“Modeling  of  Steady  Laminar  Flames  for  One-dimensional  Premixed  Jets  of 
Heptane/Air  and  Octane/ Air  Mixtures”,  (K.  G.  Harstad  and  J.  Bellan),  AIAA-201 2- 
0340,  presented  at  the  50*  Aerospace  Sciences  Meeting,  Nashville,  TN,  January  9-12, 
2012 


11 


EXECUTIVE  SUMMARY 


The  highlights  of  the  results  from  the  manuscripts  in  Appendices  1-5  are  here  summarized. 

The  challenge  of  modeling  turbulent  reactive  flows  is  so  considerable  that  the  activity  has 
traditionally  been  decomposed  into  its  two  essential  parts:  kinetics  and  turbulence.  Usually, 
modeling  of  chemical  kinetics  has  proceeded  on  a  separate  path  from  that  of  turbulence  which 
also  includes  canonical  models  for  turbulence/reaction  interaction.  The  only  constraint  to  kinetic 
modeling  was  that  it  should  be  compact  enough  to  be  computationally  efficient  when  included  in 
a  complex  turbulent  combustion  code.  However,  there  are  definite  advantages  on  approaching 
chemical  kinetic  modeling  in  a  similar  manner  to  turbulent  flow  modeling  because  if  the 
concepts  are  similar,  the  hope  is  that  the  models  will  mesh  better  and  the  results  will  be  easier  to 
understand.  This  was  the  approach  taken  in  this  study.  The  spirit  of  the  chemical  kinetic 
modeling  approach  is  that  of  Large  Eddy  Simulations  (LES)  in  which  kinematic-energy 
significant  flow  scales  are  computed  and  the  others  are  modeled;  in  turbulence,  the  large  flow 
scales  constitute  the  former  category  and  the  small  flow  scales  constitute  the  latter  category.  The 
chemical  kinetics  parallel  is  to  obtain  a  model  which  only  retains  the  thermodynamic-energy 
significant  chemical  scales  as  progress  variables,  and  models  the  influence  of  the  other  scales. 
But  the  parallel  approach  between  kinetics  and  turbulence  was  here  extended  even  further.  In 
turbulence  modeling,  a  common  methodology  is  to  assess  the  behavior  of  the  modeled  scales  by 
analyzing  databases  created  using  Direct  Numerical  Simulations  (DNS)  in  which  all  flow  scales 
are  computed;  indeed,  current  experimental  diagnostics  do  not  permit  the  same  thoroughness  of 
information  as  that  obtained  from  DNS.  The  DNS  data  is  analyzed  in  what  is  called  an  a  priori 
study  to  inquire  about  the  behavior  of  the  small  scales  and  propose  mathematical  forms  which  fit 
this  behavior.  The  a  priori  study  is  followed  by  an  a  posteriori  study  where  the  proposed 
mathematical  forms  are  inserted  into  the  model  to  evaluate  its  performance  when  compared  to 
the  DNS  database  at  the  LES  resolution. 

A  complete  analogy  between  turbulence  and  kinetics  was  here  made  by  observing  that 
kinetic  elemental  or  skeletal  mechanisms  can  serve  for  reduced  kinetics  the  role  that  DNS  serves 
for  EES,  in  which  case  reduced  kinetic  mechanisms  can  be  viewed  as  the  complement  to  EES  in 
achieving  the  goal  of  accurate  computationally-efficient  turbulent  reactive  flow  simulations.  The 
analogy  between  reduced  kinetic  models  and  EES  is  not  entirely  surprising  since  each  chemical 
species  has  a  characteristic  time  scale  and  in  the  kinetic  reduction  it  is  desirable  to  compute  only 
those  entities  (e.g.  species,  combination  of  species,  radicals,  combination  of  radicals,  etc.)  having 
essential  characteristic  time  scales  (to  be  defined)  and  model  the  kinetics  of  the  remaining 
entities. 

Thus,  there  were  two  important  components  to  this  study,  namely  the  a  priori  analysis  and 
the  a  posteriori  evaluation.  Whereas  turbulence  modeling  benefits  from  decades  of  work  using 
the  DNS/EES  concept  which  originated  in  atmospheric  turbulence  predictions  in  the  1960s,  the 
present  work  is  the  first  investigation  to  take  this  approach  in  kinetic  reduction  modeling. 
Therefore,  it  was  first  necessary  to  produce  a  categorization  of  scales  analogous  to  the  large  and 
small  scales  of  turbulence,  then  it  was  required  to  propose  mathematical  forms  for  the  scales  to 
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be  modeled  rather  than  eomputed  as  progress  variables  in  the  redueed  kinetie  model,  and  finally 
it  was  required  to  perform  an  a  posteriori  study  in  order  to  evaluate  the  ehemieal  kinetie  model 
versus  the  elemental  or  skeletal  meehanism  for  those  speeies  predieted  by  the  redueed  model. 
The  present  model  depiets  a  eonstant-volume  situation,  so  as  to  be  eonsistent  with  the 
requirement  of  a  LES  grid. 

We  have  proposed  sueh  a  eategorization  [1,2]  through  the  definition  of  a  total  eonstituent 
molar  density  whieh  is  a  progress  variable  representing  the  heavy  speeies  kineties,  and  through 
the  partition  of  the  light  speeies  set  into  a  set  of  modeled  quasi-steady  speeies  and  a  set  of 
progress  variable  speeies.  By  definition,  eonstituent  radieals  are  those  eomposing  speeies  with  a 
earbon  number  larger  than  or  equal  to  3.  These  speeies  are  ealled  ‘heavy’,  and  their  eomplement 
in  the  ensemble  of  speeies  is  ealled  ‘light’.  Constituents  are  obtained  by  mathematieally  breaking 
the  heavy  speeies  into  parts.  Although  sometimes  lights  and  eonstituents,  both  of  whieh  are 
radieals,  may  have  the  same  ehemieal  formula,  the  differenee  between  eonstituents  and  these 
light  speeies  is  that  the  later  are  unbound  to  other  ehemieal  entities  whereas  the  former  are  bound 
to  other  ehemieal  entities,  i.e.  other  eonstituents.  The  element  eompositions  of  the  eonstituents 
are  not  linearly  independent,  but  the  eonstituents  are  linearly  independent.  Thus,  the  eonstituents 
are  independent  struetural  elements  whieh  have  individual  valenee  bond  topologies;  the 
eonstituents  are  not  just  based  on  atom  eounts.  The  total  eonstituent  molar  density  is  the  sum  of 
the  individual  eonstituents’  molar  densities. 

Thus,  rather  than  following  all  speeies  through  their  reaction  coordinates,  we  follow  a 
reduced  set  of  reaction  coordinates  (i.e.  progress  variables);  this  reduced  set  is  called  a  base.  An 
extensive  analysis  of  the  LENT  databases  [3]  of  elementary  sets  of  reactions  for  a  given  fuel 
species  revealed  that  a  normalized  temperature  can  be  defined  which  serves  as  a  similarity 
variable,  0.  The  variable  which  is  self-similar  versus  0  is  the  molar  density  of  the  constituents 
divided  by  the  product  of  the  equivalence  ratio  and  a  reference  nitrogen  molar  fraction  which 
serves  as  a  surrogate  pressure.  Both  the  a  priori  and  a  posteriori  models,  with  numerous  results, 
were  described  for  n-heptane  [1],  for  iso-octane,  for  mixtures  of  iso-octane  and  n-heptane  (which 
are  called  Primary  Reference  Euels)  and  for  mixtures  of  iso-octane  and  n-pentane  or  of  iso¬ 
octane  and  iso-hexane  [2].  The  reduced  models  were  compared  with  the  original  elementary 
EENL  mechanisms  [3].  Comparisons  involved  temperature,  species  molar  densities,  the 
constituent  molar  density  and  ignition  times.  The  results  were  uniformly  excellent  except  for  the 
very  large  equivalence  ratios  (i.e.  4)  where  they  were  only  very  good.  The  deterioration  in 
performance  for  these  very  rich  flames  was  due  to  the  multivalued  aspects  of  the  functions 
versus  0. 

Specifically,  we  have  shown  that  the  molar  density  of  the  constituents  divided  by  the  product 
of  the  equivalence  ratio  and  a  reference  nitrogen  molar  fraction  plotted  versus  0  is  invariant  with 
pressure  over  the  range  5-50  atm,  with  the  equivalence  ratio  over  the  range  [1/4,  4]  range  and 
with  initial  temperature  values  in  the  cold  ignition  regime,  i.e.  [600  K,  900  K].  Eor  larger  initial 
temperatures  than  900  K  and  up  to  1200  K  to  which  the  model  has  been  tested,  the  self-similarity 
holds  modulo  the  initial  temperature  value.  Remarkably,  this  model  is  valid  over  a  much  larger 
equivalence  ratio  range  than  all  other  chemical  kinetic  reduced  mechanisms  presented  in  the 
literature  that  are  typically  developed  in  the  [1/2,  2]  range.  As  part  of  another  research  program. 
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recent  simulation  results  describing  the  mixing  of  five  species  under  supercritical  pressure  [4] 
show  that  the  [1/2,  2]  range  captures  an  extremely  limited  range  of  locations  in  the  flow  field.  An 
additional  feature  of  the  self-similarity  parameter  0  is  that  the  molar  densities  of  oxygen  and 
water  plotted  versus  0  display  a  quasi  linear  decreasing  and  quasi  linear  increasing  behavior, 
respectively,  over  the  entire  range  of  parameters  explored.  The  fact  that  these  features  prevail  for 
not  only  single-fuel  species  but  also  mixtures  of  fuel-species,  including  mixtures  of  different 
percentages  of  n-heptane  and  iso-octane,  is  very  encouraging  for  its  extension  to  heavier 
hydrocarbons.  In  fact,  our  model  is  hierarchical  by  construct,  meaning  that  it  is  naturally 
extendable  to  higher  carbon-number  hydrocarbons  (a  primary  necessity  for  modeling  heavier 
than  iso-octane  and  ring  hydrocarbons),  thus  having  advantages  in  this  respect  over  the 
capabilities  of  other  models.  The  hierarchical  aspect  comes  in  because  the  progress  variables  are 
not  necessarily  species,  but  also  include  species  ‘constituents’  defined  very  much  like  in  group 
additivity  theory  [5].  The  number  of  constituents  in  the  global  constituent  may  increase  with 
increasing  carbon  number  or  with  moving  from  straight-chain  alkanes  to  ring  hydrocarbons,  but 
the  self-similarity  is  expected  to  remain.  An  example  of  the  increasing  number  of  constituents 
has  already  been  observed  on  going  from  n-heptane  to  iso-octane  (13  versus  14),  but  the  primary 
pillar  of  the  method  which  is  the  self-similarity,  still  holds. 

To  give  perspective  to  our  model,  the  number  of  species  progress  variables  in  our  model  is 
11,  whereas  reduced  models  in  the  scientific  literature  claiming  good  agreement  with  data  have 
more  than  50  progress  variables.  The  saving  in  computational  time  using  our  model  is  very 
significant. 

On  consulting  chemical  physicists  [6]  in  retrospect,  there  seems  to  be  no  surprise  to  our 
success.  This  is  because  the  heavy  species  decompose  very  fast  upon  heating,  and  it  is  only  what 
results  from  this  decomposition  that  makes  an  impact  on  the  reaction.  We  have  taken  advantage 
of  this  property  to  develop  a  compact,  yet  accurate  chemical  kinetic  model. 

Building  on  the  success  with  chemical  kinetics  reduction,  the  next  topic  addressed  has  been 
that  of  laminar  flame  propagation.  Compared  to  chemical  kinetics  prediction  alone,  flames  have 
the  additional  complication  of  involving  transport  of  species,  momentum  and  energy.  Therefore, 
accurate  transport  coefficients  must  be  computed  for  the  mixtures  under  consideration.  For 
example,  the  species  mass  diffusion  matrix  is  a  square  matrix  having  as  many  elements  in  each 
of  the  two  directions  as  species  and  species-like  progress  variables.  Each  element  of  the  matrix 
depends  on  the  composition,  temperature  and  pressure  [7,  8].  The  thermal  diffusion  factor  matrix 
is  involved  in  the  computation  of  Soret  and  Dufour  effects  which  may  be  important  under  high- 
pressure  conditions  [8].  The  thermal  conductivity  must  be  computed  for  the  entire  mixture,  with 
appropriate  mixing  rules  accounting  for  the  varying  composition  and  temperature  during  the 
evolution  of  the  flame.  Although  the  flow  is  considered  inviscid,  according  to  [9]  the  viscosity 
plays  a  role  in  the  computation  of  thermal  conductivity,  so  it  must  be  modeled.  Additionally,  a 
real-gas  equation  of  state,  having  capability  to  perform  well  under  high-pressure  conditions, 
must  also  be  constructed,  with  appropriate  mixing  rules  [9].  Moreover,  because  of  the 
partitioning  of  the  species  into  constituents,  quasi-steady  lights  and  progress-variable  lights,  the 
governing  equations  required  now  recasting  in  a  new  form  to  respect  this  partition.  This 
preliminary  work  for  computing  flame  has  been  presented  in  [10]  and  applied  to  computing  the 
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evolution  of  the  temperature  and  speeies  in  a  premixed  flame,  up  to  the  flame  loeation,  using  a 
eomputational  method  based  on  a  large  Peclet  number.  To  assess  the  importance  of  species  mass 
diffusion,  computations  were  performed  with  and  without  species  mass  diffusion  and  the  effect 
of  species  mass  diffusion  has  been  highlighted.  Further,  by  changing  the  numerical  method  to  a 
split-step  type  where  there  are  many  chemical  steps  within  a  diffusion  numerical  step, 
computations  were  carried  out  for  premixed  flames  through  the  flame  region  [11].  Currently,  this 
work  is  being  enlarged  to  perform  simulations  for  preparing  a  manuscript  to  be  submitted  for 
publication  in  a  refereed  journal.  The  next  step  will  be  the  establishment  of  a  methodology  for 
computing  counterflow  flames  within  the  constituent  concept. 
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A  methodology  for  deriving  a  reduced  kinetic  mechanism  for  alkane  oxidation  is  described  and  applied  to 
n-heptane.  The  model  is  based  on  partitioning  the  species  of  the  skeletal  kinetic  mechanism  into  lights, 
defined  as  those  having  a  carbon  number  smaller  than  3,  and  heavies,  which  are  the  complement  in  the 
species  ensemble.  For  modeling  purposes,  the  heavy  species  are  mathematically  decomposed  into  con¬ 
stituents,  which  are  similar  but  not  identical  to  groups  in  the  group  additivity  theory.  From  analysis  of 
the  LLNL  skeletal  mechanism  in  conjunction  with  CHEMKIN  II,  it  is  shown  that  a  similarity  variable 
can  be  formed  such  that  the  appropriately  scaled  global  constituent  molar  density  exhibits  a  self-similar 
behavior  over  a  very  wide  range  of  equivalence  ratios,  initial  pressures  and  initial  temperatures  that  is  of 
interest  for  predicting  n-heptane  oxidation.  Furthermore,  the  oxygen  and  water  molar  densities  are 
shown  to  display  a  quasi-linear  behavior  with  respect  to  the  similarity  variable.  The  light  species  ensem¬ 
ble  is  partitioned  into  quasi-steady  and  unsteady  species.  The  concept  is  tested  by  using  tabular  informa¬ 
tion  from  the  LLNL  skeletal  mechanism  in  conjunction  with  CHEMKIN  11.  The  test  reveals  that  the 
similarity  concept  is  indeed  Justified  and  that  the  combustion  temperature  is  well  predicted,  but  that 
the  ignition  time  is  overpredicted.  To  palliate  this  deficiency,  functional  modeling  is  incorporated  into 
our  conceptual  reduction.  Due  to  the  reduction  process,  models  are  also  included  for  the  global  constit¬ 
uent  molar  density,  the  kinetics-induced  enthalpy  evolution  of  the  heavy  species,  the  contribution  to  the 
reaction  rate  of  the  unsteady  lights  from  the  heavies,  the  molar  density  evolution  of  oxygen  and  water, 
the  mole  fractions  of  the  quasi-steady  light  species  and  the  mean  molar  heat  capacity  of  the  heavy  spe¬ 
cies.  The  model  is  compact  in  that  there  are  only  nine  species-related  progress  variables.  Results  are  pre¬ 
sented  comparing  the  performance  of  the  model  for  predicting  the  temperature  and  species  evolution 
with  that  of  the  skeletal  mechanism.  The  model  reproduces  the  ignition  time  over  a  wide  range  of  equiv¬ 
alence  ratios,  initial  pressure  and  initial  temperature. 

©  2010  The  Combustion  Institute.  Published  by  Elsevier  Inc.  All  rights  reserved. 


1.  Introduction 

The  reduction  of  elementary  or  skeletal  oxidation  kinetics  to  a 
subgroup  of  tractable  reactions  for  inclusion  in  turbulent  combus¬ 
tion  codes  has  been  the  subject  of  numerous  studies.  The  skeletal 
mechanism  is  obtained  from  the  elementary  mechanism  by  remov¬ 
ing  from  it  reactions  which  are  considered  negligible  for  the  intent 
of  the  specific  study  considered.  As  of  now,  there  are  many  chemical 
reduction  methodologies.  A  typical  way  of  reducing  chemical 
mechanisms  is  to  represent  it  by  a  few  significant  (in  quantitatively 
defined  ways)  chemical  species  and  correspondingly  perform  a 
lumping  of  Arrhenius-type  reactions.  Examples  of  models  that  fall 
in  this  category  are  those  of  Mueller  et  al.  [1  ],  Bollig  et  al.  [2],  Sung 
et  al.  [3]  and  Li  et  al.  [4].  Another  method,  called  piecewise  imple¬ 
mentation  of  solution  mapping  (PRISM)  [5],  relies  on  the  construc¬ 
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tion  of  a  polynomial  representation  in  a  region  of  chemical 
composition  space  to  approximate  the  chemical  kinetics.  This  rep¬ 
resentation  is  stored  for  use  during  a  numerical  simulation,  and 
the  chemistry  evolution  for  points  within  the  composition  space 
is  computed  by  evaluating  these  polynomials  instead  of  solving  dif¬ 
ferential  equations.  If  the  evolution  of  chemistry  is  required  outside 
of  the  domain  where  the  polynomials  are  available,  the  methodol¬ 
ogy  dynamically  samples  the  region  and  constructs  a  new  represen¬ 
tation.  In  situ  adaptive  tabulation  (ISAT)  [6]  is  a  method  similar  to 
PRISM,  and  also  relies  on  tabulated  polynomial  fits  to  the  chemical 
mechanism  over  a  given  time  interval  at  various  locations  in  the 
composition  space.  ISAT  is  based  on  a  linear  approximation, 
whereas  PRISM  represents  a  piecewise  quadratic  fit  of  the  chemis¬ 
try.  The  intrinsic  low-dimensional  manifold  (ILDM)  [7]  method  is 
based  on  the  observation  that  in  a  chemically  reactive  mixture  at 
fixed  pressure,  fixed  total  enthalpy  and  fixed  atomic  element  com¬ 
position,  the  path  of  reaction  in  the  high-dimensional  state  space 
lies  in  low-dimensional  manifolds  in  the  vicinity  of  reaction  attrac- 
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tors,  far  from  the  final  equilibrium  point;  these  manifolds  are 
intrinsic,  meaning  that  except  for  the  degree  of  reduction,  the  reac¬ 
tion  mechanism  in  terms  of  Arrhenius-type  elementary  reactions  is 
the  only  source  of  information  necessaiy  to  build  the  approxima¬ 
tion.  However,  ILDM  is  not  easily  implementable  in  some  situa¬ 
tions.  Particularly,  it  (i)  has  the  drawback  of  becoming  very 
difficult  to  handle  with  increasing  number  of  C  atoms  in  the  fuel 
[8],  and  we  note  that  heavy  species  (i.e.  large  number  of  C  atoms) 
enter  the  composition  of  most  practical  fuels,  and  (ii)  does  not  work 
at  low  temperatures,  or  very  lean  or  rich  mixture  conditions,  which 
are  precisely  the  conditions  where  a  veiy  large  number  of  pollu¬ 
tants  are  formed.  The  Computational  Singular  Perturbation  (CSP) 
method  of  Lam  and  Goussis  [9]  introduces  simplifications  by  distin¬ 
guishing  between  the  different  time  scales  of  the  elementary  reac¬ 
tions.  The  CSP  method  computerizes  the  choice  of  the  important 
reactions,  and  therefore  does  not  rely  on  intuition  to  retain  the 
important  slow  reactions  with  respect  to  the  fast  ones.  However, 
CSP  reduces  only  the  number  of  reactions,  not  the  number  of  spe¬ 
cies  [9],  That  is,  the  number  of  differential  equations  to  be  solved 
In  order  to  find  the  evolution  of  the  mixture  remains  the  same  as 
in  the  detailed  mechanism,  and  thus,  although  rigorous  and  accu¬ 
rate,  CSP  remains  computationally  expensive.  Lumping  procedures, 
both  for  reactions  and  for  components  have  also  been  used  [10,11] 
for  mechanism  reduction;  the  components  are  usually  possible 
Isomers  of  large  hydrocarbons  and  In  this  procedure  a  large 
number  of  real  components  is  lumped  into  a  judiciously  selected 
number  of  equivalent  components.  Thus,  the  lumped  mecha¬ 
nisms  of  heavy  species  are  represented  by  a  limited  number  of 
equivalent  reactions.  The  directed  relation  graph  (DRG)  reduction 
[12-16]  is  a  more  recent  method  in  which  the  coupling  among 
species  is  plotted  In  a  graph  which  is  then  analyzed  to  identify 
unimportant  species;  these  species  are  then  removed  from  the 
mechanism.  These  few  reduction  methodologies  only  represent 
a  subset  of  all  procedures  devised  for  the  goal  of  obtaining  a  ki¬ 
netic  mechanism  compact  enough  to  be  utlllzable  for  turbulent 
reactive  flow  calculations  and  accurate  enough  to  be  reliable 
over  a  wide  range  of  equivalence  ratios,  </>,  initial  pressures,  pg 
(subscript  0  denotes  the  Initial  value),  and  initial  temperatures. 
To.  As  none  of  the  existing  methodologies  for  mechanism  reduc¬ 
tion  is  considered  the  ultimate  answer  in  the  quest  for  compact¬ 
ness  and  reliability,  the  search  for  novel  ideas  in  kinetic 
mechanism  reduction  continues.  The  study  presented  here  is 
the  result  of  such  a  search. 

In  this  study,  one  of  the  main  concerns  was  to  devise  a  kinetic 
reduction  procedure  which  is  consistent  with  the  Large  Eddy  Sim¬ 
ulation  (LES)  concept,  as  LES  has  shown  considerable  promise  In 
simulating  turbulent  flows  and  represents  the  state-of-the-art 
capability  in  such  simulations.  The  primary  idea  in  LES  is  that 
due  to  the  computational  infeasibillty  of  routinely  solving  the  gov¬ 
erning  equations  for  fully  turbulent  flows  (i.e.  high  Reynolds  num¬ 
ber),  one  should  spatially  Alter  the  equations  thus  removing  the 
dynamic-energy  unimportant  small  scales,  leaving  only  the  large 
scales  to  be  resolved.  The  Influence  of  the  small  scales  is  re-intro- 
duced  in  the  equations  through  functional  modeling,  allowing  the 
solution  of  the  large  scales,  and  thus  of  most  of  the  dynamic  en¬ 
ergy.  The  dynamic  scales  of  fluid  mechanics  have  a  parallel  In 
chemical  kinetics  since  each  species  has  its  own  characteristic  time 
and  can  be  thought  to  be  a  scale  of  the  problem.  The  dynamic  en¬ 
ergy  of  the  flow  has  a  parallel  in  the  thermodynamic  energy  of  the 
reaction;  at  this  early  stage  of  our  reduction  model,  we  are  only 
Interested  in  recovering  the  energetics  of  the  reaction,  the  major 
species,  and  some  of  the  species  through  which  a  reaction  is  mon¬ 
itored  (e.g.  OH).  Therefore,  the  idea  here  is  to  explore  whether  it  is 
possible  to  remove  scales  unimportant  to  the  energetics,  and  func¬ 
tionally  model  them.  The  study  focusses  on  the  cold  ignition  re¬ 
gime  for  hydrocarbon  oxidation. 


In  LES,  the  small  scales  are  modeled  by  examining  solutions  ob¬ 
tained  from  Direct  Numerical  Simulation  for  incipient  (i.e.  low  Rey¬ 
nolds  number)  turbulence,  as  these  solutions  contain  small-scale 
behavior  information.  For  chemical  kinetics,  the  information  that 
could  be  neglected  is  contained  in  the  skeletal  kinetic  mechanism. 
Therefore,  the  idea  is  here  that  examination  of  the  solutions  from 
the  skeletal  mechanism  should  reveal  the  functional  forms  of  the 
kinetic  scales  which  could  be  neglected. 

This  paper  is  organized  as  follows;  We  first  describe  the  ba¬ 
sis  leading  to  our  conceptual  model.  Then,  we  present  the  mod¬ 
el  which,  because  it  results  from  examination  of  the  n-heptane 
kinetic  mechanism  [17],  is  directly  intertwined  with  the  quanti¬ 
tative  aspects  of  n-heptane  oxidation.  The  interest  is  in  n-hep- 
tane  because  it  is  the  simplest  hydrocarbon  exhibiting  a 
negative  temperature  coefficient  (NTC)  behavior  and  is  a  refer¬ 
ence  fuel.  Further,  we  assess  the  model  by  comparing  it  with 
the  skeletal  mechanism  to  identify  modeling  needs.  The  func¬ 
tional  model  Is  next  presented,  and  results  from  it  are  critically 
examined  for  different  aspects  of  the  predictions.  In  a  following 
section,  we  discuss  the  ensemble  of  the  results  to  indicate  the 
strategy  for  future  work.  Finally,  a  summary  and  conclusions 
are  presented. 

2.  Conceptual  model 

The  conceptual  model  seeks  to  represent  the  species,  thought  to  be 
akin  to  vectors  in  a  mathematical  space,  by  a  reduced  base  set.  The 
species  are  partitioned  into  heavies  (carbon  number,  n  >  3)  and 
lights  (the  remaining  of  the  set).  The  heavies  can  be  either  radicals 
or  stable  species.  The  lights  are  oxygen,  nitrogen,  the  final  combustion 
products  and  light  radicals/molecules  (e.g.  CH3,  CH4,  H2O2).  The 
heavies  are  decomposed  into  constituents,  as  described  below,  while 
the  lights  remain  part  of  the  base. 

2.2.  Heavy  species:  the  total  constituent  molar  density  as  progress 
variable 

The  definition  of  constituents  stems  from  the  observations  that 
(i)  plots  of  the  heats  of  combustion  for  alkanes  and  for  alkenes  hav¬ 
ing  a  C  double  bond  at  the  molecular  chain  end  have  a  linear  var¬ 
iation  with  the  C  number,  n  and  (ii)  at  fixed  temperature  T,  the 
species  molar  heats  at  constant  pressure,  Cp,  vary  linearly  with  n. 
The  implications  are  that  (1)  heats  of  combustion  and  Cp  may  be 
considered  obtainable  by  summing  those  of  constituent  radicals 
CH2,  CH3  and  C2H3  that  form  these  hydrocarbons  and  (2)  for 
n  ^  3,  species  may  be  mathematically  decomposed  into  constitu¬ 
ent  radicals.  The  mathematical  decomposition  is  not  meant  to 
emulate  a  real  reaction,  such  as  pyrolysis.  Each  constituent  molar 
density,  Nk,  is  the  sum,  over  all  heavy  species,  of  the  count  of  the 
constituent  in  each  heavy  species  multiplied  by  the  molar  density 
of  that  species.  The  element  compositions  of  the  constituents  are 
not  linearly  Independent,  but  the  constituents  are  linearly  inde¬ 
pendent.  Basically,  the  constituents  are  independent  structural  ele¬ 
ments  which  have  valence  bond  topologies  (see  below);  the 
constituents  are  not  just  based  on  atom  counts. 

The  constituents  are  similar  to  the  groups  in  group  additivity 
theory  [18,19],  however,  there  are  marked  differences  from  it.  In 
group  additivity  one  accounts  for  interactions  with  adjacent  groups, 
interactions  with  non-adjacent  groups  and  for  steric  effects.  In  our 
‘constituents’  concept  we  only  account  for  interactions  with 
adjacent  groups  and  first  order  (compositional)  effects.  Also,  our 
process  is  different  from  lumping  because  we  are  decomposing  all 
heavy  species  and,  as  explained  above,  a  constituent  may 
span  the  entire  species  set  of  heavy  species.  The  heavy  species 
constituent  radicals  are  defined  as  a  set  of  entitles  (here,  13  of 
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Fig.  1.  Structure  of  the  13  constituents.  First  row,  left  to  right:  methylene,  methyl, 
formyl.  Second  row,  left  to  right:  vinyl,  vinylene.  Third  row,  left  to  right:  dicarbon, 
ethynyl,  keto.  Fourth  row,  left  to  right:  methylidyne,  hydroxyl,  hydroperoxy.  Fifth 
row,  left  to  right:  peroxy,  atomic  oxygen. 

them:  CH2,  CH3,  CH,  C2H3,  C2H2,  C2,  HC2,  CO  (keto),  HCO,  HO, 
HO2, 00,  0)  from  which  any  heavy  species  or  radical  may  be  con¬ 
structed,  The  structural  topology  of  the  constituents  is  illustrated 
In  Fig,  1,  This  mathematical  decomposition  is  precise  for  the  large 
majority  of  heavy  species.  If  not  precise,  a  particular  constituent 
which  is  in  the  heavy  species  but  not  in  our  set  of  13  constituents, 
is  replaced  by  the  closest  constituent  in  the  set.  The  total  constitu¬ 
ent  molar  density  is 

(1) 

li:=l 


where  the  molar  density  of  individual  constituent  k  Is  Nh. 

Each  case  computed  Is  Initialized  at  a  specified  pg,  Tg  and  tj). 
Replacing  the  pressure  variable  p,  we  define  a  dimensionless  molar 
density 


IV*  = 


ref 


(2) 


for  N2,  where  the  reference  N2  molar  density  is  Nref  =  31.5  mol/m^ 
for  dry  air  value  at  pressure  p„f  =  1  bar  and  !„/  =  298.15  K.  Thus, 
N*  acts  as  a  convenient  surrogate  variable  for  p  since  the  N*  value 
Is  essentially  reaction  rate  invariant  given  that  N2  is  inert  and  no 
NOx  mechanism  is  considered  at  this  early  stage  of  model  develop¬ 
ment.  That  Is,  the  partial  pressure  of  N2  is  the  overwhelming  contri¬ 
bution  to  p;  basically,  the  ratio  of  N*  to  pg  is  nearly  constant. 

We  seek  to  obtain  a  normalized  variable  0  such  that  at  approx¬ 
imately  the  same  value  of  0,  Nc  has  been  consumed  In  the  reaction 
for  all  initial  conditions,  except  for  rich  situations.  If  such  a  variable 
can  be  obtained,  one  would  hope  to  formulate  the  problem  in  more 
generic  terms  than  if  this  variable  could  not  be  found.  Beside  the 
goal  of  achieving  a  normalization  so  that  Nc  decreases  by  3  orders 
of  magnitude  at  approximately  same  0  (except  for  rich  situations), 
we  also  ask  the  question  whether  the  normalization  variable  could 
additionally  be  a  similarity  variable  such  that  all  Nc  versus  6  curves 
nearly  coincide.  The  advantage  of  self-similarity  is  that  it  reduces 
the  dimensionality  of  the  problem.  Because  of  the  complexity  of 
hydrocarbon  oxidation  mechanisms,  performing  a  general  mathe¬ 
matical  analysis  to  seek  6  Is  unpractical,  so  the  search  for  9  was 
performed  by  choosing  a  specific  example,  the  n-heptane  oxidation 
LLNL  skeletal  mechanism,  and  by  exhaustively  examining  the  data¬ 


base  through  curve  fitting.  The  result  of  this  very  time-consuming 
empirical  work,  described  in  Section  3.1,  was 


p-  r-To 

Tr  =  2065(N*)°“w((^) 


w(<^)  =  (j) 


1.5  +  1.31.^ 

TTcTtT^TTT.^' 


(3) 

(4) 

(5) 


As  it  will  be  shown  in  the  following,  the  9  definition  is  such  that  in¬ 
deed  Nc  decreases  by  three  orders  of  magnitude  (delving  into  higher 
order  of  magnitude  decrease  runs  the  risk  of  encountering  round-off 
and  truncation  errors)  from  its  initial  value  during  stoichiometric 
combustion  upon  reaching  9  >  0.6.  The  9  ~  0.6  value  was  chosen 
to  ensure  that  all  9  values  remain  below  unity  for  all  test  case 
calculations. 


2.2.  Light  species:  a  modeled  subset  and  a  progress  variable  subset 

According  to  the  above  procedure,  light  species  are  not  subject 
to  meaningful  decomposition.  Examination  of  runs  made  using  the 
LLNL  database  in  CHEMKIN  II  Indicates  that  the  light  species 
should  be  categorized  In  two  subsets:  one  which  Is  modeled,  and 
one  which  is  computed  subject  to  the  heavy  species  model,  as 
follows. 


2.2.1.  Quasi-steady  light  species 

The  species  in  the  first  subset  are  the  radicals  O,  CH, 
CH2,  CH3,  HO,  HCO,  HO2,  HC2,  C2H3  which  have  a  quasi-steady 
behavior  (production  and  consumption  rates  are  within  5%  of  each 
other  during  an  overwhelming  time  of  the  reaction).  Their  mole  frac¬ 
tion,  Xj,  is  here  modeled  through  mathematical  fits  as  a  function  of 
the  state  variables  (()'),N*,r)  and  modeling  parameters  (here,  only 
To).  There  are  nine  of  these  quasi-steady  light  species. 


2.2.2.  Progress  variable  light  species 

The  species  in  the  second  subset  are  unsteady  and  require  rate 
equations.  This  set  consists  of  H2O,  CO2,  O2,  H,  CO,  H2,  CH4, 
H2O2,  C2H2,  C2H4,  CH2O.  The  reaction  rate  of  light  species  is  con¬ 
ceptualized  as 


fdNQ 

_dNi 

^dNi 

\dtj 

reac 

heavies 

lights 


(6) 


where  the  first  term  expressing  the  contribution  to  the  lights  from 
of  the  heavy  group  must  be  modeled  and  second  term  In  the  right 
hand  side  has  the  same  rates  as  in  the  LLNL  skeletal  mechanism. 
We  make  the  ad-hoc  assumption  that  to  lowest  order,  through 
our  reduction 


dN, 

dt 


13 


J^NQKGg,-XiKLglc) 

k=^ 


(7) 


where  KCih  and  XL,  ^  represent  gain  and  loss  rates  from  the  heavy 
constituent  k  to  the  light  i.  Examination  of  the  skeletal  mechanism 
result  reveals  that  individual  dominant  constituents  of  molar  den¬ 
sity  Nji  are  quasi-steady,  which  implies  that  there  is  a  mole  fraction 
XCi!  such  that  Nj,  ~  NcXCk,  meaning  that 


dN, 

dt 


■Nc 


^XQKGi_,-X,^XCMi.i 


(8) 


By  defining  XG,  =  ZlijNQKGck  and  XL,  ee  Y^H^XCkKLgk,  one  obtains 


dN,' 


dt 


heavies 


Nc(KGi-XiKL,), 


(9) 


where  XG,  and  XL,  are  functions  of  (ro,<^,N*,T)  and  they  must  be 
modeled  consistently  with  the  heavy  species  model.  There  are  ele¬ 
ven  of  these  unsteady  light  species. 
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2.3.  Model  summary  for  species  and  computation  of  energetics 

At  this  junction,  the  base  set  is  composed  of  the  total  constitu¬ 
ent  radical  molar  density  Nc  and  of  the  1 1  molar  densities  of  the 
light  molecules  or  free  radicals  of  the  second  light  species  subset. 
Thus,  there  are  at  this  stage  only  12  species-related  progress  vari¬ 
ables.  To  use  this  model,  one  must  first  find  a  strategy  for  comput¬ 
ing  Nc\  such  a  strategy  will  be  shown  next.  To  compute  the  20  light 
species,  one  must  model  Xj  of  the  first  light  species  subset  and 
compute  the  conventional  light  species  reaction  rates  of  the  second 
light  species  subset.  For  the  unsteady  light  species  reaction  rates, 
models  of  XG,  and  KL,  are  needed  in  functional  form. 

Computation  of  the  temperature  evolution  in  a  reactive  system 
requires  knowledge  of  the  species  molar  enthalpies  and  molar  heat 
capacities^  For  species  i,  the  molar  enthalpy  may  be  expressed  as 
hi  =  hf  +  hi(T)  where  h°  is  the  heat  of  formation,  h,-  is  the  sensible 
enthalpy  and  T  is  the  temperature.  The  heat  of  formation  is  taken  at 
the  above  reference  conditions,  giving  hi  =  Cp  idT  where  Cp  ,  is 
the  molar  constant-pressure  heat  capacity.  Under  the  assumption 
that  water  remains  in  vapor  state,  a  heat  of  combustion  for  species 
i  is  given  by 

h'i^h°-Y,Wph°  (10) 

j 

where  Wj,-  is  the  number  of  species  j  in  species  i  that  produces  the 
correct  atom  count,  where  the  index  j  denotes  the  species  set  of 
air  {O2,  N2)  and  final  combustion  products  (H2O,  CO2).  (This  means 
that  mathematically,  we  consider  species  i  as  being  composed  of  air 
and  final  combustion  product  species.)  Values  of  fi)  and  h°  are  pro¬ 
vided  in  Table  1. 

Considering  the  total  enthalpy  to  be  fixed,  and  consistent  with 
the  goal  of  only  modeling  the  energetics  resulting  from  the  oxida¬ 
tion  reaction,  the  temperature  evolution  is  given  by 


Table  1 

Thermodynamic  properties  of  molecules  and  free  radicals.  and  h'  (heats  of 
formation  and  combustion,  respectively),  in  kj/mol;  constants  a''  and  b'’  for  molar  heat 
capacity  in  the  form  Cp/R„  =  a*'  +  b''ln(T/T„f);  T„j  =  298.15  K.  “Mo"  denotes  “mol¬ 
ecule”.  “Ra"  denotes  “radical”. 


Mo/Ra 

X 

a'’ 

b'' 

H2 

0.0 

241.5 

3.282 

0.400 

O2 

0.0 

0.0 

3.476 

0.5663 

N2 

0.0 

0.0 

3.388 

0.469 

c 

717 

nil 

2.50 

0.0 

H20 

-241.5 

0.0 

3.688 

1.217 

C02 

-393.5 

0.0 

4.690 

1.390 

N 

473 

473 

2.50 

0.0 

H 

218.0 

339 

2.50 

0.0 

HO 

38 

159 

3.385 

0.3637 

HOO 

10.5 

131 

4.150 

1.307 

0 

249.2 

249.2 

2.536 

0.0 

CO 

-110.5 

283 

3.426 

0.4749 

HCO 

43.1 

558 

4.154 

1.2875 

CH4 

-74.6 

802 

3.797 

4.305 

CH3 

146 

902 

4.440 

2.249 

CH2 

390 

1025 

3.973 

1.3015 

CH 

596 

1110 

3.220 

0.7136 

C2H3 

300 

1449 

5.1 

3.5 

HC2 

566 

1474 

4.434 

1.404 

C2 

838^ 

1625 

4.58 

0.0 

NO 

90.3 

90.3 

3.533 

0.4508 

NO2 

33.2 

33.2 

4.691 

1.151 

H2O2 

-136 

106 

5.269 

1.880 

HCOH 

-109 

526.5 

4.27 

2.546 

C2H2 

228 

1257 

5.368 

2.294 

C2H4 

52.5 

1323 

5.383 

4.676 

®  Alternate  value  of  832  from  the  CRC  tables. 


where  Ru  is  the  universal  gas  constant;  this  equation  is  general  and 
holds  for  the  skeletal  mechanism.  Equation  (11)  is  an  approxima¬ 
tion  of  the  energy  equation  in  that  it  does  not  include  the  tran- 
sient-p  term,  so  as  not  to  constrain  the  validity  of  the  results  to 
only  low-p  situations  for  which  the  perfect-gas  equation  of  state 
(EOS)  used  in  CHEMKIN  is  appropriate.  CHEMKIN  11  contains  two 
non-equivalent  options  for  calculations  of  reaction  rates:  (1)  pro¬ 
viding  T  and  the  molar  densities,  a  procedure  which  does  not  neces¬ 
sitate  the  use  of  the  EOS,  and  (2)  specifying  p,  additionally  to  T  and 
the  molar  densities,  a  procedure  which  requires  the  utilization  of 
the  perfect-gas  EOS.  Since  high-p  combustion  will  necessitate  the 
use  of  real-gas  EOS,  in  our  computations,  we  use  option  (1),  so  as 
to  make  our  results  independent  of  the  EOS.  The  neglect  of  the  p 
variation  term  in  Eq.  (11)  is  justified  because  (i)  during  the  initial 
phase  before  ignition,  the  T  variation  is  very  small  and  the  corre¬ 
sponding  changes  in  p  are  also  small,  and  (ii)  during  ignition,  when 
the  timewise  p  variation  is  significant,  most  of  the  T  increase  is  due 
to  the  heat  release  from  combustion.  Neglect  of  the  transient-p  ef¬ 
fects  corresponds  to  slow  and/or  small  changes  in  the  mean  local 
mass  density  (this  is  consistent  with  fixed  total  enthalpy  results  ob¬ 
tained  using  the  LLNL  skeletal  model  in  CHEMKIN  11;  not  shown),  i.e. 
an  approximately  incompressible  fluid. 

For  the  reduced  model,  since  the  heavy  species  are  not  available, 
the  sums  over  the  heavy  species  must  be  replaced  by  functional 
forms  which  must  be  determined  from  numerical  calculations  with 
the  LLNL  skeletal  mechanism.  To  this  end,  we  define  a  mean  heavy 
species  molar  heat  capacity 

/-  _  i^teheavies^P.l^l) 

Cp.h  = - - ,  (12) 

and  a  mean  rate  of  enthalpy  change  (in  units  s  ') 


1 

^uT'refh^c 


(13) 


In  Eqs.  (12)  and  (13),  Nc  is  used  used  for  normalization  purposes, 
which  is  consistent  with  Eq.  (9)  and  the  conceptual  model  described 
in  Section  2.1.  Replacing  Eqs.  (12)  and  (13)  in  Eq.  (11)  yields  the 
temperature  evolution  equation 


E 


ielights 


dT 

dt 


^  hp^i  +  Nc(RuTref)I<i, 
ielights 


(14) 


in  which,  when  using  the  reduced  model,  Nc,  Cpj,  and  K/,  (all  of 
which  describe  the  heavy  species  set)  are  quantities  not  directly 
computable  from  the  reduced  model,  and  thus  must  be  modeled 
as  a  function  of  (To,  i/),N*,T);  additionally,  one  must  obtain  N,-  using 
the  model  described  in  Section  2.2.2. 

For  the  energetics,  values  of  h'^  are  obtained  from  the  literature 
and  for  each  light  species  Cp  is  modeled  as 

E  ^  (si)  ■  1’=' 

Values  of  a'’  and  fa*'  for  the  light  molecules  are  given  in  Table  1  [20-23  ]. 


3.  Results 

The  results  are  here  presented  from  four  perspectives.  First,  in 
Section  3.1  we  examine  the  benefits  of  the  reduction  of  the  heavy 
species  to  Nc  and  of  the  normalization  of  the  database  using  9.  The 
kinetic  reduction  from  the  skeletal  mechanism  to  our  model  neces¬ 
sarily  involves  a  loss  of  information  which  requires  additional 
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modeling.  Second,  to  separate  the  reduction  concept  from  the  nec¬ 
essary  additional  model,  in  Section  3.2  we  evaluate  the  concept 
using  an  ideal  model  extracted  from  the  skeletal  kinetics.  Third, 
in  Section  3.3  we  present  the  a  priori  study,  which  involves  model¬ 
ing  of  the  information  lost  through  the  kinetic  reduction  by  func¬ 
tional  fits  and  comparing  them  to  the  skeletal  mechanism 
derived  functions.  Finally,  in  Section  3.4,  in  an  a  posteriori  study, 
we  show  the  performance  of  the  kinetic  reduction  using  the  func¬ 
tional  fits. 

3.1.  Examination  of  the  n-heptane  skeletal  mechanism  using  our 
concept 

Fig.  2  shows  Nc/((i>  x  N")  as  a  function  of  0;  the  plots  were 
obtained  using  the  LLNL  skeletal  kinetics  in  CHEMKIN  11.  The 
similarity  achieved  with  Nc/{4>  x  N")  and  9  is  notable  despite 
small  departures  from  the  nominal  curves.  Even  for  as  rich  a 
mixture  as  =  2,  similarity  holds,  making  this  similarity  valid 
in  realms  beyond  those  in  advanced  reduced  schemes  where 
the  upper  limit  of  the  scheme  validity  is  </>  =  1.5  [13],  at  most 
0  =  2.0  [24]  or  exceptionally  0  =  3  [25[.  At  0  =  4,  the  mixture 
is  too  rich  for  the  reaction  to  obtain  complete  fuel  burning, 
and  as  a  result,  the  reaction  termination  induces  the  disparity 
from  self-similarity  at  0  ~  0.24. 

The  self-similar  behavior  can  be  understood  as  a  reduction  in 
dimensionality  of  the  problem  but  should  not  be  confused  with 
the  ILDM  method  [7]  because  that  concept  was  developed  for 
actual  chemical  species,  whereas  our  findings  are  for  Nc.  Basically, 
Nc  serves  as  a  coarse-grained  attractor  compared  to  the  attractors 
found  through  the  ILDM.  And  neither  is  our  model  equivalent  with 
lumping  [10,11]  because  we  are  decomposing  all  heavy  species 
rather  than  a  selection  of  them,  and  because  a  constituent  may 
span  the  entire  species  set  of  heavy  species.  The  fact  that  Nc  can 
embody  the  evolution  of  all  heavy  species  is  consistent  with,  and 
can  be  traced  to,  the  fact  that  as  T  increases,  the  heavy  species 
chemically  decompose  and  no  longer  play  a  role;  instead,  the 
products  of  this  decomposition  determine  the  reaction  evolution. 

Once  values  of  0  ~  0.6  are  achieved,  the  constituents  have  been 
nearly  exhausted  and  the  reaction  may  be  considered  to  have 
reached  completion.  Thus,  the  normalization  achieved  through  0 
may  be  very  useful  because  over  the  wide  range  of  4>,Po  and  ^o. 
the  reaction  is  completed  at  approximately  the  same  0  value,  as  de¬ 
sired  (see  Section  2.1). 

Plots  of  the  molar  densities  of  oxygen,  Noi.  and  of  water,  N/,20, 
versus  0  are  illustrated  in  Fig.  3  over  the  same  wide  range  of 
0,Po  and  To  as  in  Fig.  2.  Notably,  both  types  of  plots  display  a  qua¬ 


si-linearity.  For  H2O,  eventually,  an  asymptotic  behavior  is  seen  at 
0  ~  0.6  indicative  of  the  reaction  having  reached  completion,  con¬ 
sistent  with  the  information  from  Fig.  2. 

The  above  findings  suggest  that  instead  of  computing  Nc  from  a 
rate  equation  (i.e.  as  a  progress  variable),  one  may  simply  fit  the 
Nc/(0  X  N*)  curves  shown  in  Fig.  2;  this  fitting  must  be  performed 
in  the  four  parameter  space  (To,  0,N*,r).  That  is,  whereas  the  LLNL 
skeletal  mechanism  depends  on  more  than  the  four  variables  we 
have  chosen,  the  reduced  model  only  depends  on  these  four 
parameters;  this  fact  is  shown  next  in  Section  3.2.  Noteworthy, 
Fig.  2  shows  that  the  accuracy  of  the  Nc/(0  x  N*)  fits  versus  0  for 
0  >  0.6  is  irrelevant,  since  either  Nc  ~  0  past  this  0  value  or  for  very 
rich  situations  Nc  is  small  well  before  0  ~  0.6.  Similarly,  instead  of 
considering  O2  and  IT2O  as  progress  variables,  it  seems  reasonable 
to  functionally  fit  the  slope  of  the  curves  in  Fig.  3  and  use  this  fit  in 
routine  calculations  to  predict  these  two  major  species.  Thus,  this 
examination  of  our  conceptual  model  for  n-heptane  using  the  LLNL 
skeletal  kinetics  in  CITEMKIN  11  leads  to  further  reducing  the  num¬ 
ber  of  species-related  process  variables  from  12  to  9.  These  9  pro¬ 
gress  variables  are  the  molar  densities  of  the  lights;  CO2,  CO, 
H,  H2,  CH4,  H2O2,  C2H2,  C2H4,  CH2O. 

3.2.  Assessment  of  the  concept's  predictive  capabilities  using  an  ideal 
model 

Before  launching  into  the  task  of  performing  functional  fits  for 
rates  and  X,,  we  deemed  it  instructive  to  assess  the  predictive  capa¬ 
bility  of  the  concept,  independent  of  the  functional  fits,  by  utilizing 
ideal  functional  fits  extracted  from  the  LLNL  skeletal  mechanism 
using  CHEMKIN  11.  These  ideal  functions  were  obtained  in  tabular 
form,  every  5  °K,  and  served  as  input  to  our  conceptual  model. 
Examples  of  the  results  are  portrayed  in  Fig.  4  for  both  Nc  and  T 
profiles;  indeed,  predicting  the  energetics  and  more  particularly 
the  ignition  time,  tjg„,  is  an  important  goal  of  the  model. 

Unsurprizingly  according  to  Fig.  2,  Fig.  4a  shows  that  Nc  is  excel¬ 
lently  predicted  by  our  conceptual  model.  However,  whereas  the  va¬ 
lue  of  the  combustion  T  is  also  excellently  predicted,  as  shown  in 
Fig.  4b,  the  predicted  ignition  time  is  longer  than  that  of  the  skeletal 
mechanism,  which  epitomizes  the  missing  information  in  the  re¬ 
duced  kinetic  model.  The  fact  that  the  combustion  T is  well  predicted 
for  the  largest  t  shown  in  Fig.  4b  can  be  seen  by  mentally  translating 
the  T(t)  reduced  mechanism  curves  to  match  the  corresponding 
skeletal  mechanism  ignition  time,  to  find  that  the  translated  curves 
coincide  with  those  from  the  skeletal  model.  Experimentally,  it  is  ob¬ 
served  [26]  that  tjgn  is  an  extremely  sensitive  quantity  function  of  To 
and  even  a  1  “1<  change  in  To  makes  approximately  1 0%  change  in  tign. 


Fig.  2.  Similarity  plots  of  parameter  Nc/(i/>  x  N")  versus  ff  at  (a)  Po  =  20  bar  and  (b)  (^  =  1  using  the  LLNL  skeletal  mechanism  in  conjunction  with  CHEMKIN  11.  To  is  in  K. 
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Fig.  3.  Oxygen  and  water  molar  densities  versus  9  as  extracted  from  the  LLNL  skeletal  mechanism  in  conjunction  with  CHEMKIN  11. 


Fig.  4.  Molar  density  of  constituents  (a)  and  temperature  (b)  obtained  with  our  conceptual  model  and  tabular  information  from  the  LLNL  skeletal  mechanism  using  CHEMKIN 
II.  The  symbols  are  from  the  LLNL  skeletal  mechanism  in  conjunction  with  CHEMKIN  11  (□□□  po  =  10  atm;  O  O  O  Po  =  40  atm),  and  the  lines  are  obtained  using  tabular 
information  from  the  LLNL  skeletal  mechanism  with  CHEMKIN  II  in  our  conceptual  model  ( - Po  =  10  atm;  —  Po  =  40  atm). 


Thus,  the  accurate  prediction  of  t/gn  is  a  very  Stringent  test  for  a  mod-  scales  of  the  initial  pre-ignition  process  must  be  reintroduced 

el.  From  the  modeling  viewpoint,  this  means  that  the  neglected  through  very  careful  modeling.  This  modeling  is  described  next. 
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3.3.  Functional  fits  to  emulate  the  LLNL  skeletal  mechanism 

The  functional  fits  necessary  to  complete  the  model  fall  into 
two  categories:  fits  for  the  rates  Kh,KGi,RLi  =KLi/KGi,  for  Nc,Cpj, 
and  for  the  X,  of  the  quasi-steady  light  species;  and  fits  for  the 
slopes  of  the  unsteady  quantities  N02,  and  N^o- 

3.3.  t.  Fits  for  the  rate  quantities  and  quasi-steady  light  species  molar 
fractions 

According  to  Eq.  (14),  to  recover  the  value  of  T  in  the  reduction 
scheme,  the  heavy  species  model  should  focus  on  an  appropriate 
representation  of  Cpf,  and  Kj,.  Plots  (not  shown)  of  Cp  f,  at  various 
Po  values,  calculated  using  the  LLNL  model  over  a  wide  range  of 
<l),  show  that  these  curves  exhibit  a  very  modest  variation  over 
the  range  of  strong  Nc  decay;  moderate  (i.e.  up  to  50%)  variation 
occurs  only  after  Nc  values  are  very  small  to  negligible.  This  indi¬ 
cates  that  the  T  recoveiy  is  primarily  governed  by  Ki,  as  far  as  heavy 
species  modeling  is  concerned. 

The  model  for  Kh,KGi,RLi  andX,  is  constructed  considering  these 
quantities  as  functions  of  Tor  9,  with  N*,  tj)  and  To  being  parameters. 
Quantities  K^.KGi  and  X,  are  split  into  a  veiy  low  T  region,  i.e. 
0  C  10“^,  termed  the  incubation  region  in  which  there  are  veiy  slow 
changes;  a  low-rate  portion  corresponding  to  10“^  ^  0  ^  0.2  that 
exhibits  a  maximum,  K^x,  and  a  minimum,  K^n',  and  a  high 
T,  0  >  0.2,  high-rate  portion  termed  the  fast-rate  region  in  which  a 
continuous  increase  is  commonly  observed.  Conversely,  RL,  first 
exhibits  a  minimum,  and  then  a  maximum.  Thus,  the  functional  fits 
are  performed  in  three  separate  regions,  with  connection  constraints 
between  consecutive  regions  (the  ultimate  fits  are  piecewise 
continuous). 

Curve  fits  for  9  >  10“^  are  generated  utilizing  either  one  of  the 
two  following  methods.  In  the  first  method,  one  uses  a  cubic  trans¬ 
formation,  e.g.  KGi(T)  to  y(r),  such  that  y  =  -1  at  the  maximum 
point  and  y  =  1  at  the  minimum  point.  Then  values  of  T  at  fixed 
y  values  are  fitted  in  terms  of  parameters  (To,(f>,N*);  these  values 
then  generate  the  continuous  curve  y(r)  by  interpolation  of  the 
discrete  values.  Specifically,  the  non-monotonic  10“^  ^  0  C  0.2  re¬ 
gion  behavior  is  captured  by  using  a  cubic  functional  mapping 
from  K(T)  to  y(T)  through  the  form 

K(T)  =  ^(X™  +  Kp,n)  +  ^(Kmx  -  -  3y)  (16) 

where  y(T)  =  -1  at  T  =  T^x  which  is  the  location  of  K^x  and 
y(T)  =  1  at  T  =  T^n  which  is  the  location  of  In  fact,  the  range 
of  y  is  from  slightly  smaller  than  -2  at  T  =  Q  (i.e.  temperature  at 
which  Kc  =  -d(lnNc)/dt  >  0)  to  y  <  1.8.  The  values  of 
Kmx,I(mn,T„,x  and  Tmn  ate  fitted  as  polynomials  in  parameters 
ln(N*),ln((f))  and  ln(rs/rre/).  Temperature  Ts(To,N'',(j>)  is  fitted  as  is 
T„  defined  as  the  T  value  at  which  y  =  0.  Two  different  changes  of 
variables  are  made  from  y(T)  to  y(z)  depending  on  how  T  compares 
to  T„.  For  T  <  Tn,  a  variable  z  =  (!„  -  T) /(!„  -  Tmx)  is  defined,  which 
means  that  y  =  -latz=l.To  achieve  the  mapping,  we  first  gener¬ 
ate  (y,z)  pairs  which  are  functions  of  parameters  ln(N*),ln((f))  and 
ln{rs/rre/).  To  match  these  pairs,  a  set  of  appropriate  functions 

(e.g.  polynomials  a  +  bz  +  cz^-i - ;  power  functions  of  type  a^; 

and  combinations  of  the  two  functional  forms)  is  chosen  to  produce 
the  y(z)  mapping.  For  r  >  r„,  i.e.y  >  0,  a  similar  procedure  is  used 
for  y(z)  where  now  z  =  (T  -  T„)/(Tmn  -  Tn)-  Beyond  0  5=  0.2,  for  the 
high-T  region,  ln(X)  is  fitted  as  a  fifth  order  polynomial  in  9’’  where 
p  =  1  for  (f)  ^  0.5  and  p  =  (0.86  +  O.280)/{O.74  +  O.Slf)  for  tj)  >  0.5. 
These  polynomial  coefficients  are  fitted  in  terms  of  ln{N*),  In(if))  and 
in(ro/r„^). 

In  the  second  method,  the  T  intereals  before  the  maximum 
point,  between  maximum  and  minimum  point  and  after  minimum 
point  are  treated  separately.  Each  of  these  intervals  is  divided  into 
equal  T  slices,  and  logarithm  values  at  slice  boundaries  are  func¬ 


tionally  fitted  in  terms  of  parameters  (Toji^.N*).  The  discrete  set 
of  these  equally  spaced  functional  forms  is  used  to  generate,  by 
polynomial  interpolation,  the  continuous  function.  This  second 
method  is  also  used  sometimes  for  the  high-rate  region  fits.  The 
choice  of  the  particular  method  for  any  i  is  determined  by  the  over¬ 
all  results  obtained  in  matching  the  input  functions  provided  by 
CHEMKIN  II  using  the  LLNL  database. 

The  need  to  introduce  a  model  to  capture  tig„,  as  explained  in 
Section  3.2,  as  well  as  an  extreme  sensitivity  of  /Q  on  To  (e.g. 
Kh  ~  To’  for  0  ~  10“^)  meant  that  special  care  should  be  devoted 
to  reproduce  this  dependency.  Thus,  we  developed  functional  fits 
to  the  slopes  d/Q/dT  at  To  and  made  corresponding  adjustments 
to  the  rates  at  0  <  10  ’.  These  turned  out  to  be  insufficient  to 
reproduce  tig„,  and  an  adjustment  was  introduced  to  the  initial 
slope  for  Ki,  (for  0  <  10“^)  by  using  a  multiplying  factor  determined 
by  trial-and-error  calculations,  so  as  to  best  match  the  reduced 
model  predictions  compared  to  those  of  the  skeletal  mechanism. 

This  model  was  used  for  /Q,  for  the  quasi-steady  gain  rates  XC, 
(where  i  stands  for  CO2,  H,  CO,  H2,  CH4,  H2O2,  C2H2,  C2H4, 
CH2O),  for  the  quasi-steady  loss  rate  XL,  (where  i  stands  for  H 
and  H2,  this  rate  being  null  for  the  other  light  species)  and  for 
the  quasi-steady  light  species  molar  fractions  Xi(0,  CH, 
CH2,  CH3,  HO,  HCO,  HO2,  HC2,  C2H3). 

Selected  plots  ofX/,  are  presented  in  Fig.  5;  XC^  and  XC/,  are  dis¬ 
played  in  Fig.  6  and  the  ratio  XL,  is  correspondingly  illustrated  in 
Fig.  7  for  i  =  H  and  f  =  H2.  The  rates  Xj,  exhibit  a  variation  by  as 
much  as  seven  orders  of  magnitude,  and  the  task  of  developing 
curve  fits  over  this  extended  regime  in  the  four  parameter  space 
(To,  0)  is  far  from  trivial.  Despite  the  difficulty  of  developing 
accurate  curee  fits,  the  agreement  between  the  fits  and  the  LLNL 
skeletal  mechanism  in  excellent  to  good  over  Po  6  [10,40]  bar. 
To  6  [600,800]  K  and  4>  6  [1/2,4]  and  even  extends  to  values  as 
small  as  =  1/4  (not  shown).  The  same  comments  hold  for  KG 
and  XL.  Discrepancies  between  fits  and  the  LLNL  skeletal  mecha¬ 
nism  past  0  ~  0.6  are  unimportant  for  the  reduction  concept  since 
Nc/Ncfl  <g  0(1)  past  that  station. 

For  the  quasi-steady  species,  we  chose  to  display  No  because  of 
the  high  reactivity  of  the  O  radical  and  No/,  because  the  OH  species 
is  generally  considered  as  indicative  of  the  flame  location.  The 
quantities  N,-  were  computed  as  N,-  =XjNi  where  N/,  =  fZieiightsTT’ 
thus,  the  computation  of  No  and  No/,  is  performed  in  conjunction 
with  the  calculated  total  light  species  molar  density  that  includes 
the  computation  of  the  light  unsteady  species  using  the  model  of 
Section  2.2.2  for  the  interaction  between  heavy  and  light  species. 
The  results  are  presented  in  Fig.  8  and  show  that  the  agreement  be¬ 
tween  fits  and  the  skeletal  mechanism  is  generally  very  good, 
including  some  extreme  f  values  in  the  lean  regime,  however, 
some  exceptions  occur  in  the  veiy  rich,  ^  =  A,  cases  illustrated  in 
Fig.  8e  and  f.  It  is  apparent  that  for  //>  =  4  the  model  cannot  repro¬ 
duce  the  multiple  valued  aspect  of  the  skeletal  mechanism  that  is 
due  to  the  non-monotonic  behavior  of  T(t)  for  those  cases.  We 
attribute  this  lack  of  agreement  of  the  reduced  model  with  the 
skeletal  mechanism  to  inaccuracies  in  the  fits. 

3.3.2.  Fits  for  N02  and  N/,20 

Two  quantities  that  constrain  the  evolution  of  N„2  and  N/,20  are 
the  initial  value,  Nj,2,  and  the  asymptotic  value  Nj[2o-  Thus,  fits  for 
No2  and  N/,20  were  constructed  as  follows; 

No2/Ni,2=max(1.0-A„2  X  0,0.0),  (17) 

Wfi2o/Nft2o=min(1.0,A„2o  X  0),  (18) 

and  the  slopes  A02  and  A/,20  were  fitted  as 

A„2  =  Co2(<(-,ro)x<j.x(N-)““ 

■^fj2o  =  Ch2oi4>,To)  X  {N'"f 
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Fig.  5.  Kh  versus  8  for  different  conditions.  Symbols  represent  selected  data  from  the  LLNL  runs;  lines  represent  the  corresponding  fits.  To  is  in  °1<  and  po  is  in  bar. 


where  8((p)  ^0.2  x  <j)  x  exp(-1.40)  ^  0.0525.  To  compute  functions 
Co2{<^,  To)  and  C/,2o((/>,  To),  a  set  of  (j)  values  is  first  chosen.  For  a  spec¬ 
ified  value  in  the  set,  C„2{cl)Jo)  and  Ci,2o{4>,  To)  are  fitted  as  a  polyno¬ 
mial  function  (up  to  quadratic)  of  To/Tre/.  To  obtain  values  for 
arbitrary  ij),  an  interpolation  is  performed  over  (jt  using  C02  and 
Ch2o  values  at  the  specified  0  values  in  the  chosen  set. 

Fig.  9  shows  an  example  of  the  obtained  fits  versus  the  LLNL 
skeletal  mechanism.  The  agreement  ranges  from  excellent  to 
good  and  is  typical  of  that  obtained  over  a  wide  range  of 
{To.<^.N*). 


3.4.  Model  predictions 

The  model  predictions  consist  of  the  timewise  temperature  evo¬ 
lution  r(t),  tjg„  as  extracted  from  the  T{t)  profile,  and  the  timewise 
evolution  of  the  unsteady  light  species  with  the  exception  of  O2 
and  H2O  which  are  modeled  as  shown  in  Section  3.3.2. 

3.4.1.  Temperature  profiles 

Displayed  in  Fig.  10  are  the  reduced  model  predictions  com¬ 
pared  to  those  of  the  skeletal  mechanism,  spanning  a  wide  range 
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Fig.  6.  KCi,  and  iCG^  versus  0  forpn  =  40  bar  at  different  To  and  </i.  Symbols  represent  selected  data  from  the  LLNL  skeletal  mechanism;  lines  represent  the  corresponding  fits: 
and  —  </i  =  4; <  and  —  —  tj>  =  2-,0  and - =  1;  a  and - i/i  =  l/2;>  and - <^  =  l/3;n  and  —  ■  —  c/i  =  1/4. 


of  (j)  values  and  two  values  of  each  To  and  pg-  Although  tig„  Is 
well  predicted  for  all  cases,  the  ultimate  combustion  tempera¬ 
ture  Is  not  well  predicted  at  very  lean  conditions  (<^  =  1  /4)  or 
relatively  low  pressure  (i.e.  pg  =  10  bar).  These  inaccurate  predic¬ 
tions  of  the  reduced  model  are  attributed  to  lack  of  accuracy  of 
the  mathematical  curve  fits  and  indicate,  together  with  the  con¬ 
clusions  for  0  and  OH  predictions  at  tj)  =  4,  that  perhaps  a  better 
strategy  than  fitting  would  be  to  use  the  Ideal  model  discussed 


in  Section  3.2.  As  highlighted  in  Section  3.2,  the  reduced  kinetics 
with  the  ideal  model  predicts  T  very  well;  the  comparison  with 
the  skeletal  mechanism  showed  that  only  the  ignition  time  was 
not  well  predicted.  On  the  favorable  side  for  the  results  of 
Fig.  10,  at  the  high  temperature  conditions  representative  of  die¬ 
sel,  gas  turbine  and  HCCl  engines,  the  model  performs  very  well; 
the  same  holds  for  the  very  rich  conditions  (e.g.  (ji  =  4)  at  which 
soot  forms. 
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Fig.  7.  RLft  and  RLf,2  versus  0  forpo  =  40  bar  at  different  To  and  ip.  Symbols  represent  selected  data  from  the  LLNL  skeletal  mechanism;  lines  represent  the  corresponding  fits:  0 
and  ■  ■  ■  •  =  4;fl  and  —  ■■  —  (!>  =  2;  Q  and - 0  =  1;  A  and - <l>  =  l/2;c>  and - ij>  =  \/3',D  and  —  ■  —  =  1/4. 


3.4.2.  Ignition  times 

The  success  of  the  adjusted  initial  Ki,  slopes  instrumental  in  pre¬ 
dicting  tign  is  shown  in  Fig.  11.  Noteworthy,  the  ordinate  axis  is  the 
typical  logarithmic  one  [27]  for  Fig.  11a  and  b,  but  it  is  here  linear 
for  Fig.  11c,  so  deviations  between  reduced  model  and  skeletal 
mechanism  observed  in  Fig.  11c  are  indeed  logarithmically  very 
small.  The  predictions  in  Fig.  11  are  excellent  over  a  wide  range 


of  To  including  as  low  a  temperature  as  600  K,  over  a  range  of  tf) 
including  mixtures  as  rich  as  (p  =  A,  and  as  lean  as  p  =  1/3  and 
Pa  from  1 0  to  40  bar. 

3.4.3.  Predicted  unsteady  species 

The  prediction  of  several  unsteady  light  species  is  illustrated  in 
Figs.  12  and  13.  These  results  are  obtained  from  solving  Eq.  (6)  with 
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Fig.  8.  Modeling  results  for  the  quasi-steady  molar  densities  N„  and  N„|,  and  comparison  with  those  of  the  skeletal  mechanism  at  various  To.Pg  and  ip  conditions. 


the  modeled  .3?,-.  The  Ni{6)  is  computed  from  Ni{t)  and  T{t),  which 
allows  computing  6{t). 

We  chose  to  display  N/,  and  Ni,2  because  they  both  depend  on 
the  KG  and  RL  fits,  and  thus  they  provide  a  stringent  test  of  the 
model  given  the  possibility  of  additive  error.  is  chosen  since 
CO  is  one  important  product  of  incomplete  combustion.  Ncm  is  se¬ 
lected  as  a  representative  intermediary  in  the  oxidation  process.  As 
discussed  above,  only  prediction  up  to  0  ~  0.6  should  be  consid¬ 


ered,  as  past  that  station  the  reaction  is  practically  finished  since 
Nc  ~  0  according  to  Fig.  2. 

Generally,  the  agreement  between  the  reduced  model  and  the 
skeletal  mechanism  is  good,  including  at  pg  =  10  bar  for  which 
the  combustion  temperature  was  not  well  predicted.  The  lean  re¬ 
gion  (as  lean  as  </>  =  1/3)  is  much  better  reproduced  by  the  model 
than  the  rich  region  for  which  exceptions  from  the  good  predic¬ 
tions  occur  at  (jf)  =  4  for  all  N/,,  JV/,2,  Nco  and  as  the  multivalued 
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Fig.  9.  N„2  and  N(,2„  versus  0  at  fixed  p,,  =  10  and  40  bar  for  <j>  =  l  and  To  =  700  and  800  K. 


Fig.  10.  r(t)  as  predicted  by  the  reduced  model  compared  to  the  LLNL  skeletal  mechanism  predictions  for  a  variety  of  conditions. 


region  is  not  captured  in  Fig.  12e  and  f,  or  Fig.  13e  and  f.  The  inac¬ 
curacies  at  (j)  =  4  are  conjectured  to  rise  from  the  imperfect  curve 
fits  and  indicate  again  that  the  ideal  model  instead  of  the  fits 
should  be  considered  in  the  future.  Also,  the  results  seem  to  be 
slightly  more  accurate  with  increasing  To. 

4.  Discussion 

The  goal  of  the  presented  kinetic  reduction  was  to  pragmatically 
reduce  the  skeletal  n-heptane  mechanism  to  the  smallest  possible 


set  of  progress  variables  that  still  accounts  for  the  multiple  time 
scales  occurring  in  the  NTC  regime  which  cannot  be  reproduced  by 
a  one-step  reaction.  While  current  advanced  kinetic  reduction 
schemes  having  0(50)  species-related  progress  variables  (e.g.  [28]) 
have  been  proposed  for  the  typical  range  cj)  =  [0.5, 1.0, 1.5],  compact 
reduced  n-heptane  kinetic  schemes  are  already  available  in  the  liter¬ 
ature,  but  they  also  have  only  been  illustrated  for  a  reduced  tj)  range. 
For  example,  Pitsch  and  Peters  [29[  constructed  a  reduced  n-heptane 
kinetics  scheme  for  16  species  as  progress  variables  but  for 
<!>  =  [0.5. 1.0, 2.0]  which  is  a  range  much  narrower  than  that  used 
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Fig.  11.  t|g„  as  function  ofTo,<j>  and  Pg.  Where  not  indicated  on  the  plots.  To  is  in  °I<  and  Pg  is  in  bar. 


in  the  present  study,  and  Montgomery  et  al,  [30]  utilized  in  their  re¬ 
duced  n-heptane  mechanism  from  28  to  37  progress  variable  species 
to  simulate  the  ignition  delay  for  </>  =  1 .  Our  model  only  needs  9  pro¬ 
gress  variables,  extends  to  a  much  wider  (j)  range  than  the  previous 
reduced  schemes  and  predicts  tign  very  well;  although  the  success 
in  predicting  some  species  for  very  rich  cases,  such  as  (^  =  4,  is  less 
than  desired  and  the  combustion  temperature  is  not  well  predicted 
over  the  entire  parametric  range.  This  lack  of  complete  success  is 
attributed  to  the  imperfect  mathematical  fits  over  a  very  wide  para¬ 
metric  range  involving  the  four  variables  (To,  (f),  N*,  T).  To  improve 
the  predictions,  two  avenues  are  possible  for  the  future.  The  first  pos¬ 
sibility  would  be  to  improve  the  accuracy  of  the  functional  fits.  The 
second  possibility  would  be  to  use  the  ideal  model  extracted  for 
our  concept  from  the  LLNL  mechanism  using  CHEMKIN  II.  The  incli¬ 
nation  at  this  point  would  be  to  take  the  second  avenue  both  because 
of  accuracy  and  because  the  ideal  model  would  be  very  efficient 
since  it  would  only  have  to  be  read  in  once  at  the  beginning  of  the 
computation. 

5.  Summary  and  conclusions 

A  kinetic  reduction  model  has  been  developed  that  is  based  on 
the  concept  of  constituents  representing  the  evolution  of  all  heavy 
species,  and  on  light  species  representing  the  complementary 
chemistry  to  that  of  heavy  species.  The  constituents  are  found 
through  a  mathematical  decomposition  of  the  heavy  species  and 
their  global  molar  density  is  quasi-steady.  The  fact  that  the  constit¬ 
uents  rather  than  the  heavy  species  are  important  can  be  directly 
traced  to  the  fact  that  the  heavy  species  decompose  and  it  is  the 
reaction  of  these  products  of  decomposition  that  matters  for  the 
energetics.  The  light  species  fall  into  two  categories:  quasi-steady, 
for  which  no  differential  equation  must  be  solved,  and  unsteady 


light  species  which  are  progress  variables.  A  thorough  examination 
of  the  LLNL  skeletal  mechanism  for  n-heptane  revealed  that  it  is 
possible  to  empirically  define  a  similarity  variable  which  is  func¬ 
tion  of  the  initial  temperature,  initial  pressure  and  equivalence  ra¬ 
tio,  and  for  which  a  scaled  total  constituent  molar  density  exhibits 
a  self-similar  behavior  across  initial  temperatures  in  the  cold  igni¬ 
tion  regime,  initial  pressures  and  equivalence  ratios.  The  signifi¬ 
cance  of  this  finding  is  that  there  is  a  reduction  in  the 
dimensionality  of  the  problem  and  that  the  total  constituent  molar 
density  could  be  fitted  as  a  function  of  the  representative  thermo¬ 
dynamic  variables.  Thus,  just  like  in  ILDM,  we  found  a  lower¬ 
dimensional  manifold  which  is,  however,  different  from  that  found 
in  ILDM  because  it  is  of  a  more  coarse-grained  nature.  Further 
examination  of  the  skeletal  mechanism  revealed  that  the  molar 
density  of  oxygen  and  that  of  water  displayed  a  quasi-linear  varia¬ 
tion  with  the  similarity  variable  over  the  entire  range  of  parame¬ 
ters  surveyed.  The  suggestion  was  that  the  molar  densities  of 
these  two  unsteady  light  species  could  be  functionally  modeled 
rather  than  computed  from  differential  equations.  The  equation 
determining  the  temperature  evolution  was  cast  in  the  form  corre¬ 
sponding  to  the  global-constituent  concept,  thus  identifying  the 
quantities  which  must  be  mathematically  fitted. 

The  concept  was  tested  by  using  the  LLNL  skeletal  mechanism 
in  conjunction  with  CHEMKIN  11.  That  is,  instead  of  a  model  con¬ 
sisting  of  functional  fits,  we  used  an  ideal  model  represented  by  ta¬ 
bles  extracted  from  the  LLNL  skeletal  mechanism  in  the  form 
needed  in  our  conceptual  model.  Our  conceptual  model  was  found 
to  be  sound,  but  the  ignition  time  it  predicted  was  too  long  com¬ 
pared  to  the  template.  The  meaning  of  this  discrepancy  is  that 
the  neglect  of  the  chemical  processes  during  the  initial  time  of  very 
small  temperature  changes  results  in  the  overprediction  of  the 
ignition  time.  The  situation  is  equivalent  to  turbulent  modeling 
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Fig.  12.  Predictions  of  the  reduced  model  for  Nj,  and  N/,2  compared  to  those  of  the  skeletal  mechanism  at  various  To,Po  and  0  conditions. 


using  the  Large  Eddy  Simulation  (LES)  concept  in  which  the  spatial 
filtering  of  scales  results  in  the  need  to  re-introduce  the  effect  of 
these  small  scales  through  models,  following  the  LES  modeling 
philosophy,  a  small-scale  model  was  developed  and  implemented 
at  temperatures  very  close  (1-2  °1<)  to  the  initial  temperature,  to 
be  used  with  the  conceptual  model. 

Moreover,  also  paralleling  the  LES  concept,  fits  in  functional 
form  were  developed  for  the  sensible  enthalpy  production  due  to 


the  constituents’  reaction,  for  the  quasi-steady  light  species  mole 
fractions  and  for  the  reaction  rate  of  the  unsteady  light  species. 
The  reaction  rate  of  each  unsteady  light  species  was  decomposed 
into  contributions  from  the  lights  that  were  directly  taken  from 
the  LLNL  skeletal  mechanism,  and  contributions  from  the  heavy 
species  that  were  modeled.  The  heavy  species  rate  contribution 
was  further  categorized  into  a  gain  rate  and  a  loss  rate,  which  were 
individually  modeled.  The  modeling  effort  was  by  no  means  trivial 
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Fig.  13.  Predictions  of  the  reduced  model  for  and  Nm  compared  to  those  of  the  skeletal  mechanism  at  various  To,Po  <l>  conditions. 


as  the  rates  and  molar  fractions  varied  by  as  much  as  seven  orders 
of  magnitude  and  this  variation  needed  to  be  captured  in  the  four 
parameter  space  of  the  initial  pressure,  the  initial  temperature,  the 
equivalence  ratio  and  the  temperature. 

The  consequence  of  providing  functional  fits  over  a  wide  range 
of  equivalence  ratios  encompassing  mixtures  as  rich  as  4>  =  4  and 
as  lean  as  </>  =  1  /4  was  that  generally  the  fits  were  good  to  excel¬ 


lent,  but  sparse  regions  of  only  fair  results  could  also  be  identified. 
The  ignition  times  were  excellently  to  very  well  reproduced  by  the 
model,  and  so  were  some  major  species  (e.g.  O2  and  H2  0).  Other 
major  species  (e.g.  CO)  and  OH,  which  is  an  indicator  of  the  flame 
location,  were  generally  well  reproduced  with  the  exception  of 
very  rich  situations  (i.e.  4>  =  4).  The  value  of  the  final  combustion 
temperature  was  generally  but  not  always  well  predicted.  To 
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improve  the  quality  of  the  predictions  where  lacks  were  identified, 
it  was  suggested  to  use  in  the  future  the  ideal  model  provided  by 
the  tables  extracted  from  the  LLNL  mechanism  using  CHEMKIN  11. 

Finally,  since  results  for  most  reduction  schemes  are  presented 
in  a  much  more  restricted  cj)  region  than  the  6  [1  /4, 4]  considered 
here,  it  is  difficult  to  directly  evaluate  our  model  with  respect  to 
other  models.  It  is  understood  that  for  any  given  model,  predictions 
will  improve  if  a  less  ambitious  validity  regime  is  considered. 
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A  previously  described  methodology  for  deriving  a  reduced  kinetic  mechanism  for  alkane  oxidation  and 
tested  for  n-heptane  is  here  shown  to  be  valid,  in  a  slightly  modified  version,  for  iso-octane  and  its  mix¬ 
tures  with  n-pentane,  iso-hexane  and  n-heptane.  The  model  is  still  based  on  partitioning  the  species  into 
lights,  defined  as  those  having  a  carbon  number  smaller  than  3,  and  heavies,  which  are  the  complement 
in  the  species  ensemble,  and  mathematically  decomposing  the  heavy  species  into  constituents  which  are 
radicals.  For  the  same  similarity  variable  found  from  examining  the  n-heptane  LLNL  mechanism  in  con¬ 
junction  with  CHEMKIN  II,  the  appropriately  scaled  total  constituent  molar  density  still  exhibits  a  self¬ 
similar  behavior  over  a  very  wide  range  of  equivalence  ratios,  initial  pressures  and  initial  temperatures 
in  the  cold  ignition  regime.  When  extended  to  larger  initial  temperatures  than  for  cold  ignition,  the 
self-similar  behavior  becomes  initial  temperature  dependent,  which  indicates  that  rather  than  using 
functional  fits  for  the  enthalpy  generation  due  to  the  heavy  species’  oxidation,  an  ideal  model  based 
on  tabular  information  extracted  from  the  complete  LLNL  kinetics  should  be  used  instead.  Similarly  to 
n-heptane,  the  oxygen  and  water  molar  densities  are  shown  to  display  a  quasi-linear  behavior  with 
respect  to  the  similarity  variable,  but  here  their  slope  variation  is  no  longer  fitted  and  instead,  their  rate 
equations  are  used  with  the  ideal  model  to  calculate  them.  As  in  the  original  model,  the  light  species 
ensemble  is  partitioned  into  quasi-steady  and  unsteady  species;  the  quasi-steady  light  species  mole  frac¬ 
tions  are  computed  using  the  ideal  model  and  the  unsteady  species  are  calculated  as  progress  variables 
using  rates  extracted  from  the  ideal  model.  Results  are  presented  comparing  the  performance  of  the 
model  with  that  of  the  LLNL  mechanism  using  CHEMKIN  11.  The  model  reproduces  excellently  the  tem¬ 
perature  and  species  evolution  versus  time  or  versus  the  similarity  variable,  with  the  exception  of  very 
rich  mixtures,  where  the  predictions  are  still  very  good  but  the  multivalued  aspect  of  these  functions 
at  the  end  of  oxidation  is  not  captured  in  the  reduction.  The  ignition  time  is  predicted  within  percentages 
of  the  LLNL  values  over  a  wide  range  of  equivalence  ratios,  initial  pressures  and  initial  temperatures. 

©  2010  The  Combustion  Institute.  Published  by  Elsevier  Inc.  All  rights  reserved. 


1.  Introduction 

Chemical  kinetics  simplification  has  generally  followed  several 
approaches.  Examples  are  the  approach  where  the  elementary 
mechanism  is  reduced  in  a  specified  form  which  persists  during 
an  entire  computation  (e.g.  Mueller  et  al.  [1],  Bollig  et  al.  [2],  Sung 
et  al.  [3]  and  Li  et  al.  [4]),  or  where  tabulations  storing  the  change 
in  the  state  vector  over  a  time  interval  for  later  utilization  (e.g. 
(PRISM)  [5],  (ISAT)  [6])  are  produced,  or  where  the  simplification 
in  rate  path  description  according  to  an  intrinsic  low-dimensional 
manifold  ((ILDM)  [7])  is  performed.  Kinetics  reduction  either 
through  automatic  reaction  and  species  elimination  [8]  or  through 
optimally  reduced  kinetic  models  [9]  has  also  been  proposed.  Con¬ 
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sidering  the  importance  of  accurately  predicting  the  ignition  time 
and  the  fact  that  it  is  a  very  sensitive  function  of  the  initial  temper¬ 
ature  [10],  the  goal  of  most  approaches  has  been  its  accurate  pre¬ 
diction,  although  reproducing  the  temperature  and  species 
evolution  is  also  important,  particularly  when  pollutants  are  of 
interest. 

This  study  presents  results  of  a  model  which  is  in  the  frame¬ 
work  of  both  kinetic  reduction  and  storage  of  information  for  later 
utilization,  and  proposes  a  relatively  compact  yet  accurate  model. 
This  model  is  based  on  the  kinetic  reduction  proposed  by  Harstad 
and  Bellan  [11]  through  constituents  and  species;  this  reduction 
has  been  so  far  tested  only  for  n-heptane.  Although  the  present 
model  adopts  the  same  viewpoint  and  shows  that  the  similarity 
variable  found  for  n-heptane  holds  for  iso-octane  and  its  mixtures 
with  n-pentane,  iso-hexane  or  n-heptane,  the  extension  of  the  pre¬ 
vious  model  to  higher  initial  temperatures  than  in  the  cold  ignition 
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regime  warrants  some  slight  revisions.  The  model  and  its  revisions 
are  described  in  Section  2.  Following  the  model  exposition,  is  a  pre¬ 
sentation  in  Section  3  of  the  predictions  from  the  reduced  model 
which  are  compared  for  iso-octane  and  Primaiy  Reference  Fuels 
(PRFs),  i.e.  mixtures  of  iso-octane  and  n-heptane,  to  the  corre¬ 
sponding  full  LLNL  mechanism  [12]  in  conjunction  with  CHEMKIN 
II.  Finally,  in  Section  4  a  summary  is  offered  and  conclusions  dwell 
on  the  impact  of  this  modeling  approach. 

2.  Conceptual  model 

2.1.  Summary  of  the  n-heptane  model 

In  the  conceptual  model  of  Harstad  and  Behan  [11],  the  species 
are  partitioned  into  heavies  (carbon  number,  n  >  3)  and  lights  (the 
remaining  of  the  set).  The  model  was  developed  using  the  LLNL 
skeletal  mechanism  (160  species  and  1540  reactions).  The  heavies 
can  be  either  radicals  or  stable  species.  The  lights  are  oxygen,  nitro¬ 
gen,  the  final  combustion  products  and  light  radicals/molecules 
(e.g.  CH3,  CH4,  H2O2).  The  heavies  are  not  treated  as  a  species  set, 
but  as  a  set  of  base  constituent  radicals  which  correspond  to  radi¬ 
cals  from  group  additivity  theoiy  [13[.  Although  sometimes  lights 
and  constituents,  both  of  which  are  radicals,  may  have  the  same 
chemical  formula,  the  difference  between  constituents  and  these 
light  species  is  that  the  later  are  unbound  to  other  chemical  enti¬ 
ties,  whereas  the  former  are  bound  to  other  chemical  entities,  i.e. 
other  constituents.  For  each  constituent  k,  its  molar  density,  N^, 
is  the  sum,  over  all  heavy  species,  of  the  count  of  the  constituent 
in  each  heavy  species  multiplied  by  the  molar  density  of  that  spe¬ 
cies.  The  element  compositions  of  the  constituents  are  not  linearly 
independent,  but  the  constituents  are  linearly  independent.  Thus, 
the  constituents  are  independent  structural  elements  which  have 
individual  valence  bond  topologies;  the  constituents  are  not  just 
based  on  atom  counts. 

For  each  heavy  species,  and  thus  for  the  entire  set  of  heavy  spe¬ 
cies,  there  is  a  unique  and  complete  set  of  constituents.  The  subset 
of  constituents  that  have  a  veiy  small  contribution  to  the  total  con¬ 
stituent  count  having  molar  density  Nc,  is  replaced  by  tbe  comple¬ 
mentary  subset  of  constituents  already  accounted  for  in  Nc.  Tbe 
rule  for  replacement  is  that  each  constituent  having  a  small  contri¬ 
bution  to  Nc  is  replaced  by  one  having,  first,  close  valence  structure 
and  second,  close  composition  to  the  constituent  neglected.  After 
this  replacement,  one  obtains  an  optimal  N^  set.  For  n-heptane,  the 
following  thirteen  constituents,  CH2,  CH3,  CH,  C2H3,  C2H2,  C2,  HC2, 
CO  (keto),  HCO,  HO,  HO2,  00,  0,  constitute  the  optimal  constituent 
set  and 

Nc^f^N,.  (1) 


T-Tq 

“Tr(<^,N*)’ 

Tr  =  2065(N*)°“w((j>), 


W((/>)  =  (j) 


(3) 

(4) 

(5) 


as  the  result  of  the  search,  where  T  is  the  temperature.  To  is  the  ini¬ 
tial  temperature  (subscript  0  denotes  the  initial  condition)  and  tj)  is 
the  equivalence  ratio.  It  was  shown  by  Harstad  and  Bellan  [11]  that 
indeed  Nc  decreases  by  three  orders  of  magnitude  upon  reaching 
0  >  0.6,  which  ensures  that  all  0  values  remain  below  unity  for 
all  calculations.  Moreover,  for  a  wide  range  of  pc,  Tg  in  the  range 
[600,900]  relevant  to  the  cold  ignition  regime,  and  tj)  s  [1/4,4],  all 
curves  Ncl{4>  x  N*)  as  a  function  of  0  nearly  coincided,  showing  a 
self-similar  behavior.  This  self-similar  behavior  was  exploited  in 
that  mathematical  fits  of  this  curve  replaced  the  otherwise  needed 
computations  for  NfO). 

The  light  species  were  categorized  into  a  quasi-steady  subset 
comprising  the  O,  CH,  CH2,  CH3,  HO,  HCO,  HO2,  HC2,  C2H3  radicals. 
The  radicals’  mole  fractions,  X„  were  fitted  as  functions  of  the  state 
variables  ((/),N*,T)  and  the  modeling  parameter  To.  The  complemen¬ 
tary  subset  to  the  quasi-steady  species  was  that  of  the  unsteady 
species  consisting  of  H2O,  CO2,  O2,  H,  CO,  H2,  CH4,  H2O2,  C2H2, 
C2H4,  CH2O.  Because  the  LLNL  skeletal  mechanism  showed  that 
the  mole  fractions  of  O2  and  H2O  were  quasi-linear  functions  of  0 
in  the  wide  range  of  (po,</>.To)  here  investigated,  the  approach 
was  to  functionally  fit  the  slope  of  these  two  mole  fractions  versus 
0,  so  that  the  only  remaining  species  progress  variables  were  CO2, 
H,  CO,  H2,  CH4,  H2O2,  C2  H2,  C2H4,  CH2O.  The  reaction  rates  of  un¬ 
steady  light  species  were  conceptualized  as 


dNi 

dt 


dNi 

dt 


lights 


(6) 


where  the  first  term  expressing  the  contribution  to  the  lights  from 
the  heavy  group  was  modeled  and  the  second  term  in  the  right 
hand  side  used  the  same  rates  as  in  the  LLNL  skeletal  mechanism. 
Further,  for  the  first  term. 


dNi 


dt 


heavies 


Nc{KGi-XiKLi), 


(7) 


where  KGi  and  XL,  are  functions  of  (ro,(/),N*,r)  that  were  modeled 
consistently  with  the  heavy  species  model. 

Finally,  the  temperature  evolution  was  computed  from 


-EC, 

ielights 


pM  1  ^ 


hid^i  -F  Nc{RuTref)Kh, 

ielights 


(8) 


in  which  is  the  universal  gas  constant  and 


Furthermore,  the  pressure  variable  p,  was  replaced  by  a  dimen¬ 
sionless  molar  density 


IV*  = 


ref 


(2) 


for  N2,  where  the  reference  N2  molar  density  is  Nre/=  31.5  mol/m^ 
for  dry  air  value  at  pressure  Pre/=  1  bar  and  Tref=  298.15  K.  Indeed, 
for  oxidation  in  air,  N2  is  the  dominant  species  and  thus  N*  is  a  sur¬ 
rogate  variable  for  p.  The  N'  value  is  reaction-rate  invariant  since  no 
NOx  mechanism  is  considered  at  this  early  stage  of  model  develop¬ 
ment  and  N2  is  inert. 

With  the  goal  of  obtaining  a  normalized  variable  0  such  that  at 
approximately  the  same  value  of  9,  Nc  has  been  consumed  in  tbe 
reaction  for  all  initial  conditions  (except  for  rich  situations),  the 
LLNL  skeletal  mechanism  was  exhaustively  examined  to  propose 


(9) 

(10) 


were  mathematically  modeled  as  a  function  of  (To,(/),N*,T).  For  each 
light  species  Cp  was  modeled  as 


^  =  a''  +  b*'  In 


(11) 


and  values  for  a^  and  were  provided  in  [11[. 

Thus,  the  model  required  functional  fits  for  the  molar  density 
Nc,  molar  density  slopes  for  O2  and  H2O,  X,  of  the  quasi-steady  light 
species,  XG,  and  XL,  =  XL,/XG,  for  the  unsteady  light  species,  and  X/,, 
all  as  functions  of  0  in  the  parametric  space  (ro,</),N',r). 
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2.2.  Application  of  the  n-heptane  conceptual  model  to  iso-octane  and 
its  mixtures 

When  applying  the  above  described  conceptual  model  to  an¬ 
other  alkane  than  n-heptane,  the  first  question  which  must  be 
asked  is  whether  the  number  and  list  of  constituents,  particularly 
the  optimal  set  of  constituents,  remains  the  same  as  for  n-heptane. 
Examination  of  the  LLNL  full  mechanism  for  iso-octane  (857  spe¬ 
cies  and  3606  reactions)  revealed  that  there  is  an  additional  con¬ 
stituent,  the  carbon  atom  C,  which  must  be  included  in  the 
optimal  set  because  its  molar  density  was  found  to  be  significant, 
so  that  now 

Nc^f^Nk-  (12) 

li:=l 

We  show  in  Sections  2.2.1  and  2.2.3  that  by  adopting  the  defi¬ 
nition  of  Eq.  (12),  we  can  accommodate  not  only  iso-octane,  but 
also  mixtures  of  iso-octane  with  n-pentane,  iso-hexane,  and  n-hep- 
tane,  that  0  as  defined  by  Eq.  (3)  remains  a  similarity  variable  but 
now  for  Nc  as  defined  by  Eq.  (12),  and  that  the  molar  densities  of  O2 
and  H2O  also  display  a  quasi-linear  variation  with  9  (for  <!>  =  2  and  4 
the  N/,2o  quasi-linearity  only  holds  up  to  its  maximum  value),  just 
like  for  n-heptane. 

The  model  of  Section  2.1  was  tested  in  [11]  for  n-heptane  in  the 
To  range  relevant  to  cold  ignition,  but  it  is  obviously  desirable  to 
extend  this  model,  if  possible,  to  the  higher  value  To  regime.  In  Sec¬ 
tions  2.2.2  and  2.2.3,  we  show  that  an  additional,  slight  modifica¬ 
tion  to  the  model  of  Section  2.1  makes  the  modified  conceptual 
representation  also  valid  in  the  higher  than  cold  ignition  Tq  regime. 

2.2.  J.  Iso-octane  in  the  cold  ignition  regime 

In  Figs.  1  and  2,  Nditj)  x  N*),  molar  densities  of  O2  and  H2O  (N02 
and  N/,2o)  are  illustrated  from  the  LLNL  full  mechanism  [12]  in  con¬ 
junction  with  CHEMKIN  II  for  the  same  conditions  as  in  the  respec¬ 
tive  Figs.  2  and  3  of  Harstad  and  Bellan  [11]  obtained  for  n-heptane. 
Clearly,  Ndicj)  x  N*)  remains  self-similar  with  9,  and  N02  and  N/,20 
remain  quasi-linear  with  0  (for  N/,20  up  to  its  maximum  value; 
for  (j)  =  2  and  4  after  the  N/,20  maximum  value,  the  curve  is  multi¬ 
valued  with  respect  to  0).  Visual  examination  of  curves  in  Figs.  1 
and  2  compared  to  Figs.  2  and  3  of  [11]  indicates  that  the  curees 
from  the  respective  figures  are  superimposable,  confirming  the 
similarity  concept.  Considering  iso-octane  and  n-heptane,  when 


moving  on  the  Ndid  x  N*)  curves  to  increasing  0  values,  the  impli¬ 
cation  is  that  heavy  species  pyrolysis,  decomposition  and  corre¬ 
sponding  light  species  production  occur  in  similar  manner  for  the 
two  fuels. 

Further  examination  of  the  LLNL  full  mechanism  [12]  in  con¬ 
junction  with  CHEMKIN  II  revealed  that  the  molar  density  produc¬ 
tion  rate  of  O2,  H2O,  H  and  H2  from  the  heavy  species,  as  shown  by 
Eq.  (7),  is  now  governed  by  KG,-  alone,  unlike  the  finding  for  n-hep- 
tane.  Since  O2  and  H2O  are  conceptually  computable  using  the  qua¬ 
si-linear  form  shown  in  Fig.  2  and  since  in  absence  of  diffusive 
mixing  the  molar  density  of  H  is  quasi-steady,  the  n-heptane 
form  of  the  reduced  model  is  simplified  for  iso-octane  by  using 
=  (KGi  -  XjKLj)  which  are  functions  of  (To,<f>,N\T).  All  other 
unsteady  light  species  also  have  KL,  ~  0,  as  for  n-heptane. 

2.2.2.  Iso-octane  in  the  high-temperature  regime 

To  reveal  the  applicability  of  the  conceptual  model  to  higher 
than  cold  ignition  Tq  values,  results  are  displayed  in  Fig.  3a  show¬ 
ing  Nd((l>  X  N*)  versus  0  for  To  up  to  1200  K,  as  extracted  from  the 
LLNL  full  mechanism  [12]  using  CHEMKIN  11.  Clearly,  for 
To  =  1000  K  and  1200  K,  the  values  depart  from  the  narrow  band 
in  the  ToS  [600,800]  K  regime.  In  particular,  we  observe  an  in¬ 
crease  in  the  initial  slope  of  Nd((/>  x  N*)  versus  0  with  increasing 
To,  consistent  with  the  increasing  pyrolysis  rate,  as  found  by  Yu 
et  al.  [14[.  However,  similarity  still  exists  at  fixed  Tq  as  shown  in 
Fig.  3b,  where  it  is  shown  that  for  Tq  =  1200  K  the  Nd((f>  x  N*)  val¬ 
ues  plotted  versus  0  at  various  po  and  tf)  nearly  collapse  on  each 
other.  Thus,  the  similarity  concept  still  holds  at  higher  To  but  in¬ 
stead  of  obtaining  a  single  NJjc^  x  IST)  curve  versus  0,  there  is 
now  a  family  of  curves  parametrized  by  Tq.  Extension  of  the  n-hep- 
tane  model  [11]  to  incorporate  this  new  aspect  is  discussed  in 
Section  3. 

2.2.3.  Primary  reference  fuels 

A  sample  of  results  similar  to  those  shown  in  Section  2.2.1  is  pre¬ 
sented  here  to  highlight  the  fact  that  the  similarity  concept  remains 
valid  for  PRFs.  Two  mixtures  are  chosen  to  illustrate  the  applicabil¬ 
ity  of  the  concept:  (90%  iso-octane)/(10%  n-heptane)  and  (50%  iso- 
octane)/(50%  n-heptane).  Information  extracted  from  the  LLNL  full 
mechanism  [12]  using  CHEMKIN  II  is  exhibited  for  Nd(4>  x  N") 
versus  0  in  Fig.  4;  the  quasi-linear  behavior  of  N02  and  N/,20  versus 
0  still  holds  but  is  not  illustrated  for  the  purpose  of  brevity.  The  re¬ 
sults  show  that  whether  for  (90%  iso-octane)/(10%  n-heptane)  or 


Fig.  1.  Similarity  plots  of  parameter  IVc/(0  x  N  )  versus  B  for  iso-octane  oxidation  in  the  cold  ignition  regime  at  (a)  To  =  800  K.  po  =  40  bar  and  (b)  =  1  (To  is  in  K  and  po  is  in 

bar)  using  the  LLNL  full  mechanism  in  conjunction  with  CHEMKIN  11. 
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Fig.  2.  Oxygen  and  water  molar  densities  versus  S  for  iso-octane  oxidation  in  the  cold  ignition  regime  as  extracted  from  the  LLNL  full  mechanism  in  conjunction  with 
CHEMKIN  II.  To  is  in  K  and  po  is  in  bar. 


Fig.  3.  Similarity  plots  of  parameter  Ncl(d>  x  N  )  versus  0  for  iso-octane  oxidation  at  higher  than  cold  ignition  To:  (a)  po  =  40  bar,  th  =  \  and  (b)  To  =  1200 1<  using  the  LLNL  full 


mechanism  in  conjunction  with  CHEMKIN  II.  To  is  in  K  and  po  is  in  bar. 

(50%  iso-octane)/(50%  n-heptane),  the  self-similarity  holds  very 
well  over  the  entire  range  of  0,  except  for  very  rich  mixtures  for 
which  as  0  increases  there  is  a  departure  from  the  typical  behavior 
since  the  reaction  is  incomplete;  this  departure  occurs  at  smaller  0 


with  increased  tj).  Moreover,  in  Fig.  5  is  an  example  of  self-similarity 
at  higher  To  for  the  (90%  iso-octane)/(10%  n-heptane)  mixture, 
showing  that  the  behavior  identified  for  iso-octane  holds  well 
for  PRFs. 
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Fig.  4.  Similarity  plots  of  parameter  Nc/(0  x  N  )  versus  0  at  To  =  800 1<  and  po  =  40  bar  for  a  range  of  i/i  values:  (a)  for  (90%  iso-octane)/(10%  n-heptane)  oxidation  and  (b)  for 
(50%  iso-octane)/(50%  n-heptane)  oxidation.  Results  extracted  from  the  LLNL  full  mechanism  in  conjunction  with  CHEMKIN  11. 


Fig.  5.  Example  of  similarity  plot  for  NcH<l>  x  N  )  versus  0  for  (90%  iso-octane)/(10% 
n-heptane)  at  To  =  1200 1<  for  various  po  (in  bar)  and  (/>  values.  The  results  are 
obtained  using  the  LLNL  full  mechanism  in  conjunction  with  CHEMKIN  11. 


2.2.4.  Mixtures  of  iso-octane  and  n-pentane  or  iso-hexane 

Mixtures  of  iso-octane  and  n-pentane  or  iso-hexane  display  the 
same  behavior  as  iso-octane  and  PRFs.  In  the  interest  of  brevity,  the 
corresponding  plots  are  not  shown,  but  the  predictive  capability  of 
the  model  for  such  mixtures  is  addressed  in  Section  3.3. 


3.  Results 

The  extension  of  the  conceptual  model  to  higher  To  values  than 
for  cold  ignition,  showing  now  a  families  of  similarity  cun/es 
according  to  To,  leads  to  making  the  choice  between  either  corre¬ 
spondingly  extending  the  functional  fits  of  [11]  to  encompass  this 
new  To  regime  or  using  the  ideal  model  provided  by  the  LLNL  full 
mechanism  in  the  form  of  tabular  input  for  the  conceptual  model. 
We  chose  the  second  option  both  as  a  shortcut  to  investing  addi¬ 
tional  time  in  functional  fitting  and  as  a  necessarily  more  accurate 
alternative  to  the  fits.  The  ideal  model  is  read  only  once  in  the  code, 
at  the  beginning  of  the  computation,  and  is  in  tabular  form;  during 
the  computation,  instead  of  accessing  a  function  as  was  done  in 
[11],  one  accesses  a  table.  The  functions  needed  for  the  reduced 


model  are  calculated  by  cubic  interpolation  in  T  from  the  tables, 
each  obtained  at  fixed  (To,(^,N*). 

The  listed  tabular  information  consists  of  numbers  for  T,  N^, 
quasi-steady  light  mole  fractions  Xj,  Cp  /,  as  defined  by  Eq.  (9),  K„ec.i 
and  Xh,  all  computed  from  the  LLNL  full  mechanism  [12]  in  con¬ 
junction  with  CHEMKIN  11.  The  tables  are  obtained  for  fixed  values 
of  (To,(/>,N*,r)  with  specified  increments  in  T.  Particular  attention 
was  devoted  to  T  value  increments  in  the  T  ^  To  +  5  K  regime  since 
T  may  decrease  a  few  K  (especially  for  rich  mixtures)  during  the 
initial  time  when  heavy  species  pyrolysis  occurs,  and  since  the 
reaction  spends  a  relatively  long  time  in  this  T  regime.  The  first  en¬ 
try  for  any  table  is  a  t-averaged  value  over  an  interval  up  to  the  t 
value  corresponding  to  0  ~  2  x  10  ®.  Further  table  entries  in  this 
T  C  To  +  5  I<  regime  involve  two  nested  criteria  to  determine  the  T 
frequency  of  the  entries.  In  a  first  step,  table  values  were  generated 
for  T  with  the  goal  of  computing  a  table  entry  when  (T  -  To)  is  close 
to  0.1,  0.2,  0.3,  0.5, 1.0,  and  2.5  K.  For  each  output  T,  as  one  marches 
in  t  during  the  full  model  calculations,  the  corresponding  0  value 
was  calculated  and  the  second  criterion  was  based  on  the  differ¬ 
ence  between  the  0  value  corresponding  to  this  T  and  the  value 
of  9  corresponding  to  the  previous  table  value  of  T.  At  each  occur¬ 
rence  when  this  difference  had  a  value  equal  or  larger  than  10“"*,  an 
additional  table  entry  was  generated  for  (T-  To)  ^  5  K.  However, 
once  (T-  To)  ~  2.5  K  was  reached,  only  the  criterion  based  on  the 
difference  between  9  values  was  used  until  (T-To)~5I<  was 
reached.  After  (T-  To)  ~  5  K,  tables  were  produced  at  5  K  incre¬ 
ments.  The  most  frequent  table  entries  in  the  T  C  Tq  +  5  K  regime 
were  required  for  either  larger  To  values  or  lean  mixtures. 

The  reduced  model  computations  use  the  LLNL  specified  ther¬ 
modynamic  properties  for  the  light  species.  The  species  progress 
variables  are  the  eleven  unsteady  light  species  H2O,  CO2,  O2,  H, 
CO,  H2,  CH4,  H2O2,  C2H2,  C2H4,  CH2O. 

Comparisons  are  presented  below  between  the  reduced  model 
predictions  and  those  of  the  LLNL  full  mechanism  [12]  in  conjunc¬ 
tion  with  CHEMKIN  11.  The  full  model  was  exercised  on  a  single 
processor  of  a  supercomputer  which  had  a  speed  100  times  faster 
than  that  of  a  personal  computer  on  which  the  reduced  model  was 
run.  Despite  the  speed  advantage,  the  full  model  still  required  run 
times  which  were  one  to  two  orders  of  magnitude  larger  than  the 
reduced  model. 

Results  are  presented  here  from  three  perspectives:  the  recov¬ 
ery  by  the  reduced  model  of  the  molar  density  of  important 
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Fig.  6.  Molar  density  of  (a  and  b)  H,  (c  and  d)  0,  (e  and  f)  OH  and  (g  and  h)  CO  versus  9  for  iso-octane  oxidation.  (a,c,e,g)  To  =  800  K,  po  =  40  bar  at  various  0  values  and 
(b,d,f.h)  To  =  800  K,  0  =  1  at  twopo  values.  Symbols  represent  selected  data  from  the  LLNL  runs;  lines  represent  the  corresponding  reduced  model  predictions.  To  is  in  K  and  po 
is  in  bar. 


2190 


K.  Harstad.j.  Bellan /  Combustion  and  Flame  157  (2010)  2184-2197 


Fig.  7.  T{t)  for  iso-octane  oxidation  as  predicted  by  the  reduced  model  (lines)  compared  to  the  LLNL  full  mechanism  using  CHEMKIN  II  (symbols)  for  a  variety  of  conditions.  To 
is  in  K  and  po  is  in  bar. 
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Fig.  8.  tign  for  iso-octane  oxidation  as  function  of  To,  po  and  0.  Where  not  indicated  on  the  plots.  To  is  in  I<  and  po  is  in  bar. 


species,  the  reproduction  of  T(t),  and  the  prediction  of  the  ignition 
time, 

3.1.  Iso-octane 

In  Fig.  6  are  displayed  the  molar  densities  of  H,  0,  OH  and  CO  as 
a  representative  sample  of  the  results.  H  and  0  are  highly  reactive 
radicals,  OH  is  generally  used  to  indicate  the  location  of  the  flame, 
and  CO  is  a  major  product  of  incomplete  combustion.  The  results 


are  generally  excellent  over  the  wide  range  of  tjj  =  [1/4,4]  (Fig.  6a, 
c,  e  and  g)  and  the  model  even  captures  the  turnover  of  the  curves 
at  large  0  when  the  reaction  is  nearly  completed.  In  the  rich  range, 
the  model  also  captures  the  turnover  of  the  curves  for  the  incom¬ 
plete  reaction,  but  any  multivalued  aspect,  if  it  occurs,  is  not  repro¬ 
duced.  Such  lack  of  complete  fidelity  is  unavoidable  in  reduced 
models,  as  information  is  discarded  in  the  reduction  process.  The 
reduced  model  is  also  very  successful  over  the  range  of  po  investi¬ 
gated,  as  can  be  seen  in  Fig.  6b,  d,  f  and  h. 
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Fig.  9.  Recovery  by  the  reduced  mechanism  (lines)  for  (90%  iso-octane/10%  n-heptane)  of:  (a)  Nditj)  x  N  /  versus  ffforpo  =  40  bar,  (/>  =  1  and  (b)Nc  versus  8  for  To  =  700  K,  (^  =  1. 
Comparisons  are  made  with  the  LLNL  full  mechanism  in  conjunction  with  CHEMKIN  11  (symbols).  Tq  is  in  1<  and  po  is  in  bar. 


The  prediction  of  T(t)  is  illustrated  in  Fig.  7  for  the  same  condi¬ 
tions  as  those  in  Fig.  6.  The  temperature  evolution  and  the  maxi¬ 
mum  achieved  are  very  well  predicted  by  the  reduced  model 
over  the  entire  range  of  tj)  values  (Fig.  7a).  A  slight  exception  is 
the  T  drop  in  the  rich  regime  after  the  maximum  T  achieved,  which 
is  directly  connected  to  the  uncaptured  multivalued  aspect  with 
respect  to  6  discussed  in  conjunction  with  Fig.  6.  The  decrease  of 
T  at  the  end  of  rich  combustion  is  due  to  thermal  energy  loss 
resulting  from  continued  heavy  species  reactions  which  are  not 
captured  by  the  reduction  concept;  since  this  occurs  at  the  end 
of  combustion,  it  is  expected  that  prediction  of  tig„  values  will 
not  be  affected  by  the  disparity.  For  the  leanest  case  studied,  the 
final  T  is  slightly  overpredicted  and  tjgn  is  slightly  underpredicted. 

Finally,  the  tjgn  prediction  is  shown  in  Fig.  8  for  a  wide  range  of 
To,  Po  and  (ji.  In  all  cases,  the  quantitative  aspect  of  the  prediction 
ranges  from  very  good  to  excellent,  indicating  that  the  slight 
underpredictions  observed  in  Fig.  7  are  indeed  minor.  In  particular, 
the  Negative  Temperature  Coefficient  (NTC)  region  is  excellently 
captured,  and  the  ratio  of  the  predicted  value  by  the  reduced  mod¬ 
el  to  that  from  the  full  model  is  within  percentages  rather  than 
multiplying  factors. 

3.2.  Iso-octane/n-heptane  mixtures 

To  highlight  the  flexibility  of  the  model  in  computing  kinetics 
for  PRFs,  we  present  here  detailed  results  for  two  PRF  mixtures 
of  very  different  compositions. 

3.2.1.  Ninety  percent  iso-octane/ten  percent  n-heptane 

Figure  9  displays  comparisons  between  the  reduced  and  full 
models  over  a  wide  range  of  Tq  and  po  values.  For  all  runs  the  pre¬ 
dictions  of  the  reduced  and  full  models  coincide.  Of  note  is  the  Nd 
(</)  X  N*)  versus  0  behavior  according  to  families  of  Tq  (Fig.  9a)  and 
the  lack  of  self-similarity  of  alone  versus  0  (Fig.  9b);  indeed,  only 
the  normalization  NcKtj)  x  N*)  leads  to  self-similarity  versus  0.  To 
highlight  the  capabilities  of  the  reduced  model,  illustrated  in 
Fig.  10  are  the  molar  densities  of  O2  and  H2O  over  a  wide  paramet¬ 
ric  range.  In  all  cases  with  the  exception  of  the  very  rich  cases,  N02 
and  N/,2o  are  excellently  predicted.  Even  for  the  very  rich  cases,  the 
prediction  is  still  very  good  and  the  model  captures  the  turnover  in 
the  cuiri/es,  however,  the  multivalued  aspect  is  missed,  just  as  for 
iso-octane,  a  fact  which  is  again  traced  (see  below)  to  the  incapa¬ 
bility  of  the  kinetic  reduction  to  reproduce  the  decreased  T  at  the 


end  of  combustion.  The  N„  results  shown  in  Fig.  1 1  emulate  those 
for  No2  and  N/,20  in  that  the  results  are  excellent  with  the  exception 
of  the  very  rich  cases  for  which  they  are  still  very  good.  Finally,  a 
sample  of  the  No  and  Not,  predictions  is  exhibited  in  Fig.  12  showing 
again  the  high  fidelity  of  the  predictions. 

In  Fig.  13a  and  c  the  temporal  variation  of  T  is  shown  from  our 
reduced  model  for  the  (90%  iso-octane)/(10%  n-heptane)  mixture 
oxidation  and  compared  to  the  full  model  results.  Similarly  to 
the  iso-octane  results,  the  ultimate  value  of  T  is  excellently  pre¬ 
dicted,  but  there  is  an  underprediction  of  tjgn  that  increases  with 
decreasing  po;  that  is  to  say  that  the  best  predictions  are  at  high 
pressure  which  is  the  regime  of  choice  for  the  operation  of  diesel, 
gas  turbine  and  HCCl  engines.  As  stated  above,  for  the  very  rich 
cases,  the  temperature  decrease  at  the  end  of  combustion  is  not 
captured,  which  explains  the  lack  of  complete  accuracy  in  predict¬ 
ing  the  molar  densities  in  that  regime.  Finally,  tjgn  is  illustrated  for 
this  PRF  mixture  in  Fig.  14a,  c  and  e  over  a  range  of  Tq  6  [600,1200] 
1<  (the  plots  are  versus  1000/To  to  emulate  similar  plots  for  other 
reduced  kinetics  in  the  literature),  po  6  [10,50]  bar  and  (/)  =  [1/ 
4,4[.  The  general  observations  are  that  the  predictions  as  a  function 
of  To  are  excellent  and  capture  with  high  fidelity  the  NTC  region, 
and  those  versus  po  and  tj)  are  excellent  for  high  po  and  stoichiom¬ 
etric  and  lean  mixtures.  A  slight  deterioration  is  observed  for  small 
Po  and  large  tj),  but  even  in  those  cases  the  ratio  of  the  prediction  by 
the  reduced  model  to  that  by  the  full  model  is  only  within  percent¬ 
ages.  Moreover,  the  results  versus  po  are  plotted  using  a  regular 
rather  than  logarithmic  axis,  and  thus  indeed  the  deviation  from 
the  full  model  predictions  is  veiy  small. 

3.2.2.  Fifty  percent  iso-octane/fifty  percent  n-heptane 

The  results  for  the  (50%  iso-octane)/(50%  n-heptane)  mixture 
are  similar  to  those  for  the  (90%  iso-octane )/( 10%  n-heptane)  mix¬ 
ture,  and  thus  we  only  show  some  examples  to  document  the  re¬ 
duced  model  capabilities.  Presented  in  Fig.  15  are  N02  and  N/,20: 
displayed  In  Fig.  16  are  the  molar  densities  of  the  two  radicals  0 
and  H  and  the  light  species  H2;  exhibited  in  Fig.  17  are  the  molar 
densities  for  the  flame  indicator  radical  OH  and  the  product  of 
incomplete  combustion  CO.  All  evaluations  for  the  (90%  iso-oc- 
tane)/(10%  n-heptane)  mixture  remain  valid  for  the  (50%  iso-oc- 
tane)/(50%  n-heptane)  mixture. 

For  comparison  with  the  (90%  iso-octane)/(10%  n-heptane)  mix¬ 
ture,  T(t)  is  plotted  in  Fig.  13b  and  d  and  tjgn  in  Fig.  14b,  d  and  f. 
Compared  to  the  predictions  for  (90%  iso-octane )/( 10%  n-heptane). 
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Fig.  10.  Comparisons  between  reduced  mechanism  predictions  (lines)  for  (90%  iso-octane/10%  n-heptane)  and  corresponding  LLNL  full  mechanism  using  CHEMKIN  II 
(symbols)  for  the  molar  densities  of  (a.c.e)  O2  and  (b,  d,  f)  H2O  for  a  variety  of  conditions.  Figure  (a)  has  the  same  legend  as  figure  (b).  To  is  in  I<  and  po  is  in  bar. 


those  for  (50%  iso-octane)/(50%  n-heptane)  are  somewhat  less 
favorable  in  the  cold  ignition  regime  and  also  versus  the  entire 
range  of  po,  although  they  are  still  very  good;  the  same  quality  as 
for  the  (90%  iso-octane)/(10%  n-heptane)  mixture  is  maintained 
versus  0. 

As  a  last  test,  tig„  is  plotted  in  Fig.  14g  versus  the  mole  fraction  of 
iso-octane  in  the  fuel,  Xoa.  on  a  regular  (rather  than  logarithmic) 
axis  and  the  predictions  become  better  with  decreasing  Xoa, 
although  even  for  Xoa  =  1  they  are  still  very  good. 


3.3.  Iso-octane/n-pentane  and  iso-octane/iso-hexane  mixtures 

Since  species  and  temperature  predictions  are  less  sensitive  to 
modeling  errors  than  ignition  times,  for  brevity,  we  only  focus  here 
on  tjgn  since  any  inaccuracies  in  the  model  are  magnified  when  pre¬ 
dicting  this  quantity.  In  Table  1  are  listed,  as  an  example,  ignition 
times  for  iso-octane/n-pentane  and  iso-octane/iso-hexane  mix¬ 
tures  obtained  with  the  reduced  model  and  with  the  LLNL  full 
mechanism  with  CHEMKIN  II;  for  comparison,  the  corresponding 
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Fig.  11.  Predictions  of  the  reduced  mechanism  (lines)  compared  to  the  LLNL  full  mechanism  using  CHEMKIN  II  (symbols)  for  the  (90%  iso-octane/10%  n-heptane)  mixture 

showing  the  molar  density  of  Nc©  versus  9  at  various  conditions.  The  legend  for  figure  a  is:  y  and  —  —  0  =  4;  0  and - ^  =  2:0  and  —  (/>  =  !;□  and  —  •  —  (/>  =  1/2;  a  and  -  - 

-  0  =  1  /4.  To  is  in  K  and  po  is  in  bar. 


Fig.  12.  Predictions  of  the  molar  density  for  radicals  O  and  OH  by  the  reduced  mechanism  (lines)  for  (90%  iso-octane/10%  n-heptane)  compared  to  the  corresponding  LLNL  full 
mechanism  in  conjunction  with  CHEMKIN  II  (symbols).  Selected  conditions.  To  is  in  K.  po  is  in  bar. 


iso-octane  case  is  listed  as  well.  Clearly,  the  predictions  of  the  re-  to  those  of  the  full  model.  This  performance  is  typical  of  the  suite  of 
duced  model  are  within  approximately  20%  or  better  with  respect  results. 
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Fig.  13.  for  PRF  fuels:  (a  and  c)  (90%  iso-octane)/(10%  n-heptane)  and  (b  and  d)  (50%  iso-octane)/(50%  n-heptane)  at  various  conditions.  Predictions  of  the  reduced  model 
compared  those  of  the  LLNL  mechanism  used  with  CHEMKIN  11.  Figure  (a)  has  the  same  legend  as  figure  (b).  Tq  is  in  K  and  po  is  in  bar. 


4.  Summary  and  conclusions 

A  model  previously  developed  for  reducing  n-heptane  oxidation 
kinetics  has  been  here  slightly  modified  to  address  the  oxidation  of 
iso-octane,  Primary  Reference  Fuels  (PRFs)  and  iso-octane  mixtures 
with  n-pentane  and  iso-hexane.  The  conceptual  model  was  based 
on  the  identification  of  constituents  as  the  building  blocks  of  the 
heavy  species  and  on  reducing  the  species  progress  variables  to  fol¬ 
lowing  the  total  constituent  molar  density  and  the  molar  densities 
of  the  unsteady  light  species.  A  normalized  surrogate  temperature 
variable  had  been  defined  for  which  a  scaled  total  constituent  mo¬ 
lar  density  exhibited  a  self-similar  behavior  for  initial  tempera¬ 
tures  in  the  cold  ignition  regime  and  for  a  wide  range  of  initial 
pressures  and  equivalence  ratios;  also  the  molar  densities  of  O2 
and  H2O  exhibited  a  quasi-linear  behavior  versus  the  normalized 
variable.  For  n-heptane,  mathematical  fits  versus  the  normalized 
surrogate  temperature  variable  were  required  in  a  four-dimen¬ 
sional  space  to  solve  the  equations  and  obtain  the  species  and  tem¬ 
perature  evolution.  On  going  from  the  previous  study  of  n-heptane 
to  the  present  study  of  iso-octane  and  its  mixtures,  examination  of 
the  LLNL  full  mechanisms  using  CHEMKIN  II  showed  that  the  con¬ 
stituent  list  was  augmented  by  one  additional  entity,  but  the  molar 
density  of  this  total  constituent  remained  self-similar  versus  the 
previously  defined  normalized  surrogate  temperature  variable  for 
initial  temperature  values  in  the  cold  ignition  regime,  and  for  the 
wide  range  of  initial  pressures  and  equivalence  ratios  previously 
examined.  When  the  higher  than  cold  ignition  temperature  regime 
was  explored,  a  family  of  self-similar  curves  versus  the  normalized 


surrogate  temperature  variable  was  obtained  according  to  the  ini¬ 
tial  temperature.  The  molar  densities  of  O2  and  H2O  remained  qua¬ 
si-linear  versus  the  normalized  surrogate  temperature  variable. 

Based  on  these  observations,  it  was  concluded  that  the  effort  of 
extending  the  mathematical  fits  developed  for  n-heptane  was  not 
warranted,  and  instead  an  ideal  model  in  tabular  form  was  used, 
as  extracted  from  the  LLNL  full  mechanisms  using  CHEMKIN  II 
and  our  concept.  This  tabular  information  is  read  into  the  code  only 
once,  at  the  beginning  of  the  computation,  and  is  accessed  simi¬ 
larly  to  the  mathematical  fits  for  n-heptane.  There  are  eleven  spe¬ 
cies  progress  variables,  as  the  molar  densities  of  O2  and  H2O  which 
were  obtained  for  n-heptane  using  fits  for  their  slopes,  are  now 
computed  in  an  unsteady  way  utilizing  the  ideal  model. 

Results  were  presented  from  computations  using  the  reduced 
model  for  the  oxidation  of  iso-octane,  of  two  different  PRF  fuels, 
and  of  mixtures  of  iso-octane  and  n-pentane  or  iso-hexane.  These 
results  were  compared  to  the  corresponding  information  extracted 
from  the  LLNL  full  mechanisms  using  CHEMKIN  11.  The  compari¬ 
sons  were  performed  for  initial  temperatures  in  the  [600,1200]  K 
range,  initial  pressures  in  the  [5,50]  bar  regime  and  equivalence  ra¬ 
tios  encompassing  mixtures  with  air  as  rich  as  <^  =  4  and  as  lean  as 
4i  =  'i-14.  The  capability  of  the  model  was  assessed  from  the  view¬ 
points  of  predicting  species,  the  temperature  evolution  and  the 
ignition  time.  The  species  prediction  was  excellent  in  all  cases, 
with  the  exception  of  the  multivalued  regions  occurring  for  very 
rich  situations.  The  temperature  values  were  correspondingly 
excellently  to  very  well  predicted,  except  for  the  temperature  de¬ 
crease  from  the  maximum  observed  for  rich  mixtures.  Finally, 
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Fig.  14.  tjgn  for  PRF  fuels:  (a,c,e)  (90%  iso-octane)/(10%  n-heptane)  and  (b.d.f)  (50%  iso-octane)/(50%  n-heptane)  versus  1000/To,  Po  and^;  (g)  versus  the  mole  fraction  of  iso¬ 
octane  in  the  fuel,  Xoct-  Where  not  indicated  on  the  plots,  To  is  in  I<  and  po  is  in  bar.  Comparisons  are  made  between  predictions  of  the  reduced  model  and  those  of  the  LLNL 
mechanism  used  with  CHEMKIN  11. 
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Fig.  15.  Molar  density  of  (a)  O2  and  (b)  H2O  for  (50%  iso-octane)/(50%  n-heptane)  obtained  from  the  reduced  model  (lines)  and  the  LLNL  mechanism  using  CHEMKIN  11 
(symbols):  v  ^nd  —  •  ■  —  0  =  4;  ^  and - <f>  =  2\0  and  —  0  =  1;  n  and  —  ■  —  0  =  1/2;  A  and - 0  =  1/4.  To  is  in  1<  and  po  is  in  bar. 


Fig.  16.  Predictions  of  the  molar  density  of  (a)  0,  (b)  H  and  (c)  H2  for  (50%  iso-octane)/(50%  n-heptane)  from  the  reduced  model  (lines)  and  the  LLNL  mechanism  used  with 
CHEMKIN  II  (symbols)  for  selected  conditions.  To  is  in  I<  and  po  is  in  bar. 


the  ignition  times  were  reproduced  within  percentages  of  those 
computed  with  the  full  mechanism. 

The  favorable  predictive  aspect  of  our  model  should  be  con¬ 
trasted  with  its  relative  simplicity,  as  it  has  a  specified  form  which 


is  not  dynamically  changed  during  the  computation  and  only  re¬ 
quires  solutions  of  eleven  species  progress  variables,  all  of  which 
are  light  species;  the  model  has  much  increased  range  of  </>  validity 
with  respect  to  other  reduced  models  with  a  small  number  of 
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Fig.  17.  Predictions  of  the  molar  density  of  (a  and  c)  OH  and  (b  and  d)  CO  for  (50%  iso-octane/50%  n-heptane)  by  the  reduced  model  (lines)  and  the  corresponding  LLNL  full 
mechanism  used  with  CHEMKIN  11  (symbols)  at  various  conditions.  To  is  in  K  and  po  is  in  bar. 


Table  1 

Ignition  times  for  oxidation  of  iso-octane  and  iso-octane/n-pentane  or  iso-octane/iso- 
hexane  mixtures.  Comparisons  between  the  LLNL  full  mechanism  used  with 
CHEMKIN  11  and  our  reduced  model.  To  =  800  K,  po  =  40  bar  and  <l>  =  \. 


Fuel  composition 

LLNL 

Reduced 

tjgn  ratio 

model  tign 

(s) 

(S) 

/reduced'! 

\  LLNL  / 

Iso-octane 

7.22  X  10-^ 

6.13  X  10-3 

0.85 

95%  iso-octane/5%  n-pentane 

6.60  X  10-^ 

5.56  X  10-3 

0.84 

90%  iso-octane/10%  n-pentane 

6.00  X  10-^ 

5.04  X  10-3 

0.84 

50%  iso-octane/50%  n-pentane 

2.73  X  10-3 

2.19  X  10-3 

0.80 

20%  iso-octane/80%  n-pentane 

1.80  X  10-3 

1.40  X  10-3 

0.78 

90%  iso-octane/10%  iso-hexane 

5.80  X  10-3 

4.88  X  10-3 

0.84 

50%  iso-octane/50%  iso-hexane 

3.19  X  10-3 

2.69  X  10-3 

0.84 

species  [15,16].  The  advantage  of  such  a  simple  model  becomes 
Increasingly  significant  with  increasing  carbon  atoms  of  the  fuel 
because  of  the  proliferation  of  relatively  heavy  species  with 
increasing  carbon  number;  for  example,  the  skeletal  mechanism 
for  n-heptane  was  based  on  160  species  and  1540  reactions, 
whereas  the  full  Iso-octane  mechanism  has  857  species  and  3606 
reactions.  If  the  model  can  be  extended  to  even  heavier  species, 
such  as  n-nonane  (n-CgH2o),  n-decane  (n-CioH22),  n-dodecane  (n- 
C12H26).  n-tetradecane  (n-Ci4H3o)  and  n-hexadecane  (n-C]6H34),  it 
would  be  veiy  advantageous  for  performing  computations  for 
heavy  fuels. 
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A  methodology  for  deriving  a  reduced  kinetic  mechanism  for  alkane  oxidation  is  de¬ 
scribed,  inspired  by  n-heptane  oxidation.  The  model  is  based  on  partitioning  the  species 
of  the  skeletal  kinetic  mechanism  into  lights,  defined  as  those  having  a  carbon  number 
smaller  than  3,  and  heavies,  which  are  the  complement  of  the  species  ensemble.  For  mod¬ 
eling  purposes,  the  heavy  species  are  mathematically  decomposed  into  constituents,  which 
are  similar  but  not  identical  to  groups  in  the  group  additivity  theory.  Prom  analysis  of 
the  n-heptane  LLNL  skeletal  mechanism  in  conjunction  with  CHEMKIN  II,  it  is  shown 
that  a  similarity  variable  can  be  formed  such  that  the  appropriately  non-dimensionalized 
global  constituent  molar  density  exhibits  a  self-similar  behavior  over  a  very  wide  range  of 
equivalence  ratios,  initial  pressures  and  initial  temperatures  that  is  of  interest  for  predict¬ 
ing  n-heptane  oxidation.  Furthermore,  the  oxygen  and  water  molar  densities  are  shown  to 
display  a  quasi-linear  behavior  with  respect  to  the  similarity  variable.  The  light  species  en¬ 
semble  is  partitioned  into  quasi-steady  and  unsteady  species.  The  reduced  model  is  based 
on  concepts  consistent  with  those  of  Large  Eddy  Simulation  in  which  functional  forms 
are  used  to  replace  the  small  scales  eliminated  through  filtering  of  the  governing  equations; 
these  small  scales  are  unimportant  as  far  as  dynamic  energy  is  concerned.  Here,  we  remove 
the  scales  deemed  unimportant  for  recovering  the  thermodynamic  energy.  The  concept  is 
tested  by  using  tabular  information  from  the  n-heptane  LLNL  skeletal  mechanism  in  con¬ 
junction  with  CHEMKIN  H  utilized  as  surrogate  ideal  functions  replacing  the  necessary 
functional  forms.  The  test  reveals  that  the  similarity  concept  is  indeed  justified  and  that 
the  combustion  temperature  is  well  predicted,  but  that  the  ignition  time  is  overpredicted, 
which  is  traced  to  neglecting  a  detailed  description  of  the  processes  leading  to  the  heavies 
chemical  decomposition.  To  palliate  this  deficiency,  functional  modeling  is  incorporated 
into  our  conceptual  reduction.  This  functional  modeling  includes  the  global  constituent 
molar  density,  the  enthalpy  evolution  of  the  heavies,  the  contribution  to  the  reaction  rate 
of  the  unsteady  lights  from  other  light  species  and  from  the  heavies,  the  molar  density 
evolution  of  oxygen  and  water,  and  the  mole  fractions  of  the  quasi-steady  light  species. 
The  model  is  compact  in  that  there  are  only  nine  species-related  progress  variables.  Re¬ 
sults  are  presented  showing  the  performance  of  the  model  for  predicting  the  temperature 
and  species  evolution  for  n-heptane.  The  model  reproduces  the  ignition  time  over  a  wide 
range  of  equivalence  ratios,  initial  pressure  and  initial  temperature.  Preliminary  results 
for  iso-octane  using  the  full  mechanism  are  also  presented,  showing  encouragingly  that 
the  concept  may  be  generalized  to  other  alkanes.  The  utility  of  the  model  and  possible 
improvements  are  discussed. 


I.  Introduction 

The  reduction  of  elementary  or  skeletal  oxidation  kinetics  to  a  subgroup  of  tractable  reactions  for  inclu¬ 
sion  in  turbulent  combustion  codes  has  been  the  subject  of  numerous  studies.  The  skeletal  mechanism  is 
obtained  from  the  elementary  mechanism  by  removing  from  it  reactions  which  are  considered  negligible  for 
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the  intent  of  the  specific  stndy  considered.  As  of  now,  there  are  many  chemical  reduction  methodologies,  a 
few  exaples  being:  reduction  to  a  few  significant  reactions  (Mueller  et  al.,^  Bollig  et  al.,^  Sung  et  al.^  and 
Li  et  al"');  piecewise  implementation  of  solution  mapping  (PRISM) in-situ  adaptive  tabulation  (IS AT);® 
intrinsic  low-dimensional  manifold  (ILDM);^  the  Computational  Singular  Perturbation  (CSP)  method  of 
Lam  and  Goussis;®  lumping^®di  directed  relation  graph  (DRG)  reduction. These  few  reduction 
methodologies  only  represent  a  subset  of  all  procedures  devised  for  the  goal  of  obtaining  a  kinetic  mechanism 
compact  enough  to  be  utilizable  for  turbulent  reactive  flow  calculations  and  accurate  enough  to  be  reliable 
over  a  wide  range  of  equivalence  ratios,  (j),  initial  pressures,  po  (subscript  0  denotes  the  initial  value),  and 
initial  temperatures,  Tq.  As  none  of  the  existing  methodologies  for  mechanism  reduction  is  considered  the 
ultimate  answer  in  the  quest  for  compactness  and  reliability,  the  search  for  novel  ideas  in  kinetic  mechanism 
reduction  continues.  The  study  presented  here  is  the  result  of  such  a  search. 

In  this  study,  one  of  the  main  concerns  was  to  devise  a  kinetic  reduction  procedure  which  is  consistent  with 
the  Large  Eddy  Simulation  (LES)  concept,  as  LES  has  shown  considerable  promise  in  simulating  turbulent 
flows  and  represents  the  state-of-the-art  capability  in  such  simulations.  The  primary  idea  in  LES  is  that 
due  to  the  computational  infeasibility  of  routinely  solving  the  governing  equations  for  fully  turbulent  flows 
(i.e.  high  Reynolds  number),  one  should  spatially  filter  the  equations  thus  removing  the  dynamic-energy- 
unimportant  small  scales,  leaving  only  the  large  scales  to  be  resolved.  The  influence  of  the  small  scales  is 
re-introduced  in  the  equations  through  functional  modeling,  allowing  the  solution  of  the  large  scales,  and 
thus  of  most  of  the  dynamic  energy.  The  dynamic  scales  of  fluid  mechanics  have  a  parallel  in  chemical 
kinetics  since  each  species  has  its  own  characteristic  time  and  can  be  thought  to  be  a  scale  of  the  problem. 
The  dynamic  energy  of  the  flow  has  a  parallel  in  the  thermodynamic  energy  of  the  reaction;  at  this  early 
stage  of  our  reduction  model,  we  are  only  interested  in  recovering  the  energetics  of  the  reaction,  the  major 
species,  and  some  of  the  species  through  which  a  reaction  is  monitored  (e.g.  OH).  Therefore,  the  idea  here 
is  to  remove  scales  unimportant  to  the  energetics,  and  functionally  model  them. 

As  in  LES,  there  are  two  important  components  to  this  study,  namely  the  a  priori  analysis  and  the  a 
posteriori  evaluation.  In  the  a  priori  study  for  turbulence  modeling,  one  examines  the  behavior  of  the  samll 
scales  and  proposes  models  to  duplicate  it.  For  chemical  kinetics,  the  information  that  could  be  neglected 
is  contained  in  the  skeletal  kinetic  mechanism.  Therefore,  the  idea  is  here  that  examination  of  the  skeletal 
mechanism  should  reveal  the  functional  forms  of  the  kinetic  scales  which  could  be  neglected.  In  the  a 
posteriori  study,  one  includes  the  small-scales  models  in  LES  and  assesses  their  behavior  when  interacting 
with  the  large  scales  to  replicate  the  flow  which  is  known  by  other  means,  either  through  Direct  Numerical 
Simulations  or  from  experiments.  The  same  a  posteriori  methodology  is  used  here.  The  present  model  is  for 
a  constant-volume  situation,  so  as  to  be  consistent  with  the  requirement  of  a  LES  grid. 

This  paper  is  organized  as  follows:  We  first  describe  the  basis  leading  to  our  conceptual  model.  Then, 
we  present  the  model  which  results  from  examination  of  the  n- heptane  kinetic  mechanism,  as  n-heptane 
is  the  simplest  hydrocarbon  exhibiting  a  negative  temperature  coefficient  (NTG)  behavior  and  is  a  reference 
fuel.  Further,  we  assess  the  model  by  comparing  it  with  the  skeletal  mechanism  and  identify  modeling  needs. 
The  functional  model  is  next  presented,  and  results  from  it  are  critically  examined.  Preliminary  results  for 
iso-octane  based  on  the  full  mechanism  are  also  presented,  indicating  that  the  conceptual  model  is  extendable 
to  higher  alkanes.  Finally,  we  discuss  future  work. 

II.  Conceptual  model  for  n-heptane 

The  conceptual  model  seeks  to  represent  the  species,  thought  to  be  akin  to  vectors  in  a  mathematical 
space,  by  a  reduced  base  set.  The  species  are  partitioned  into  heavies  (carbon  number,  n  ^  3)  and  lights 
(the  remaining  of  the  set).  The  heavies  can  be  either  radicals  or  stable  species.  The  lights  are  oxygen, 
nitrogen,  the  final  combustion  products  and  light  radicals/molecules  (e.g.  GH3,  GH4,  H2O2).  The  heavies 
are  decomposed  into  constituents,  as  described  below,  while  the  lights  remain  part  of  the  base. 

A.  Heavy  species:  the  total  constituent  model  as  progress  variable 

The  definition  of  constituents  stems  from  the  observations  that  (i)  plots  of  the  heats  of  combustion  for  alkanes 
and  for  alkenes  having  a  G  double  bond  at  the  molecular  chain  end  have  a  linear  variation  with  the  G  number, 
n  and  (ii)  at  fixed  temperature  T,  the  species  specific  heats  at  constant  pressure,  Cp,  vary  linearly  with  n. 
The  implications  are  that  (1)  heats  of  combustion  and  Cp  may  be  considered  obtainable  by  summing  those 
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of  constituent  radicals  CH2,  CH3  and  C2H3  that  form  these  hydrocarbons  and  (2)  for  n  ^  3,  species  may  be 
mathematically  decomposed  into  constituent  radicals.  This  is  consistent  with  group  additivity  theory, 
however,  there  are  marked  differences  from  it.  In  group  additivity  one  accounts  for  interactions  with  adjacent 
groups,  interactions  with  non-adjacent  groups  and  for  steric  effects.  In  our  ‘constituents’  concept  we  only 
account  for  interactions  with  adjacent  groups  and  first  order  (compositional)  effects.  Also,  our  process  is 
different  from  lumping  because  we  are  decomposing  all  heavy  species  and  a  constituent  may  span  the  entire 
species  set  of  heavy  species.  The  heavy  species  constituent  radicals  are  defined  as  a  set  of  entities  (here,  13 
of  them:  CH2,  CH3,  CH,  C2H3,  C2H2,  C2,  HC2,  CO  (keto),  HCO,  HO,  HO2,  00,  0)  from  which  any  heavy 
species  or  radical  may  be  constructed.  The  total  constituent  molar  density  is 

13 

=  (1) 

k=l 


where  the  molar  density  of  individual  constituent  k  is  fV^,. 

We  define  a  reference  molar  density  for  nitrogen,  N2,  from  the  dry  air  value  at  pressure  Pref  =  1  bar  and 
Tref=  298.15  K,  giving  Nref=  31.5  mol/m®.  Using  We/j  a  dimensionless  N2  molar  density  is  defined  as 


N*  =  ^ 

Nref 


(2) 


which  acts  as  a  convenient  surrogate  variable  for  the  pressure  p  since  the  N*  value  is  essentially  rate  invariant. 
That  is,  the  partial  pressure  of  N2  is  the  overwhelming  contribution  to  p;  basically,  the  ratio  of  N*  to  po  is 
nearly  constant.  A  normalized  temperature  variable  6  is  defined  by 


T-n 

Tr{<^,N*) 


Tr  =  2065(iV*)°°®t(;((/)) 
1.5  +  1.310 


u;(0)  =  0 


1  +  0.710+1.10 


2  • 


(3) 

(4) 

(5) 


The  9  definition  is  such  that  A),  decreases  by  three  orders  of  magnitude  (delving  into  higher  order  of  mag¬ 
nitude  decrease  runs  the  risk  of  encountering  round-off  and  truncation  errors)  from  its  initial  value  during 
stoichiometric  combustion  upon  reaching  9  >  0.6.  The  0.6  value  was  chosen  to  ensure  that  all  9  values 
remain  below  unity  for  all  test  case  calculations.  The  relationships  Eq.  3-5  have  been  empirically  derived 
upon  examining  the  n-heptane  LLNL  skeletal  mechanism,  as  described  in  section  A. 


B.  Light  species:  a  modeled  subset  and  a  progress  variable  subset 

According  to  the  above  procedure,  light  species  are  not  subject  to  meaningful  decomposition.  Examination 
of  runs  made  using  the  LLNL  database  in  CHEMKIN  II  indicates  that  the  light  species  should  be  categorized 
in  two  subsets:  one  which  is  modeled,  and  one  which  is  computed  subject  to  the  heavy  species  model,  as 
follows: 


1.  Quasi-steady  light  species 

The  species  in  the  first  subset  are  the  radicals  0,  CH,  CH2,  CH3,  HO,  HCO,  HO2,  HC2,  C2H3  which  have 
a  quasi-steady  behavior  and  require  mathematical  fits  of  their  mole  fraction,  as  a  function  of  the  state 
variables  (0,  X*,r)  and  modeling  parameters  (here,  only  Tq).  There  are  nine  of  these  quasi-steady  light 
species. 


2.  Progress  variable  light  species 

The  species  in  the  second  subset  are  unsteady  and  require  rate  equations.  This  set  consists  of  H2O,  CO2, 
O2,  H,  CO,  H2,  CH4,  H2O2,  C2H2,  C2H4,  CH2O.  The  reaction  rate  of  light  species  is 


TZi  = 


\  _ 

V  dt  , 

Keac  dt 

heavies 


dN, 


dt 


lights 


(6) 
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where  the  contribution  of  the  heavy  group  to  the  rate  of  change  of  a  light  species  molar  density  Ni  is  expressed 
in  terms  of  quasi-steady  gain  and  loss  rates  KGi  and  KLi  as 


dt 


=  N,{KGi  -  X,KU). 

heavies 


(7) 


Since  KGi  and  KLi  are  functions  of  (Tq,  (/>,  N*,T),  they  must  be  modeled  consistently  with  the  heavy  species 
model.  There  are  eleven  of  these  unsteady  light  species. 


C.  Model  summary  for  species  and  computation  of  energetics 

At  this  junction,  the  base  set  is  composed  of  the  global  constituent  radical  mole  fraction  Nc  and  of  the 
11  light  molecules  or  free  radicals  of  the  second  light  species  subset.  Thus,  there  are  at  this  stage  only  12 
species- related  progress  variables.  To  use  this  model,  one  must  first  find  a  strategy  for  computing  Nf.;  such 
a  strategy  will  be  shown  next.  To  compute  the  20  light  species,  one  must  model  Xi  of  the  first  light  species 
subset  and  compute  the  conventional  light  species  reaction  rates  of  the  second  light  species  subset.  For  the 
unsteady  light  species  reaction  rates,  models  of  KGi  and  KLi  are  needed  in  functional  form. 

The  energy  equation  is 


dT 


c^p,h 


^  =  -  E  h,TZ,+N,{R^Tref)Ki, 


i^lights 


(8) 


iGlights 


where  hi  is  the  molar  enthalpy,  is  the  molar  constant-pressure  heat  capacity,  Ru  is  the  universal  gas 
constant, 

_  ^"^-^iGheavies 


Gp^h  = 


N, 


and  the  rate  of  enthalpy  ehange  due  to  heavy  species  I  chemical  kinetic  changes  is  defined  as 

Kl^heavies 


Kh  = 


RuTrefNc 


(9) 


(10) 


Clearly,  since  K^  is  associated  with  the  heavy  species,  it  must  be  modeled  as  a  function  of  (Tq,  0,  N*,T). 
For  each  light  species  Gp  is  modeled  as 

^=a>^  +  bHn(^).  (11) 

V  ref  J 

Values  of  and  for  the  light  molecules  are  given  in  Table  1,20-23  ^j^ej-e  values  for  the  the  heat  of 
formation,  /i?,  and  the  heat  of  combustion  for  species  i,  /i?,  are  also  given.  Both  and  are  needed  to 
compute  the  enthalpy. 


III.  Results  for  n-heptane 

A.  Examination  of  the  n-heptane  skeletal  mechanism  using  our  concept 

Figure  1  shows  Nc/{4>  x  N*)  as  a  function  of  9;  the  plots  were  obtained  using  the  LLNL  skeletal  kinetics 
in  CHEMKIN  II.  The  similarity  achieved  with  the  two  normalized  variables  is  remarkable  despite  small 
departures  from  the  nominal  curves.  Even  for  as  rich  a  mixture  as  (j)  =  2,  similarity  holds,  making  this 
similarity  valid  in  realms  beyond  those  in  advanced  reduced  schemes  where  the  upper  limit  of  the  scheme 
validity  is  (/>  =  1.5  (Lu  and  Law^^),  at  most  (j>  =  2.0  (Peters  et  al.2'^)  or  exceptionally  (j>  —  3  (Babushok 
and  Tsang2®).  At  ^  =  4,  the  mixture  is  too  rich  for  the  reaction  to  proceed  to  its  ultimate  conclusion,  and 
as  a  result,  the  reaction  termination  induces  the  disparity  from  self-similarity  at  0  ~  0.24.  This  self-similar 
behavior  can  be  understood  as  a  reduction  in  dimensionality  of  the  problem  but  should  not  be  confused  with 
the  intrinsic  low-dimensional  manifold^  because  that  concept  was  developed  for  actual  chemical  species, 
whereas  our  findings  are  for  Nc-  And,  as  stated  above,  neither  is  our  model  equivalent  with  lumping^*^d2 
because  we  are  decomposing  all  heavy  species  and  because  a  constituent  may  span  the  entire  species  set  of 
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heavy  species.  The  fact  that  Nc  can  embody  the  evolution  of  all  heavy  species  is  consistent  with,  and  can 
be  traced  to,  the  fact  that  as  T  increases,  the  heavy  species  chemically  decompose  and  no  longer  play  a  role; 
instead,  the  products  of  this  decomposition  determine  the  reaction  evolution.  Once  values  of  0  ~  0.6  are 
achieved,  the  constituents  have  been  nearly  exhausted  and  the  reaction  may  be  considered  to  have  reached 
fruition.  Thus,  the  normalization  achieved  through  9  may  be  very  useful  because  over  the  wide  range  of  </>, 
Po  and  To,  the  reaction  is  completed  at  approximately  the  same  9  value. 

Plots  of  the  molar  densities  of  oxygen,  No2,  and  of  water,  Nh2o,  versus  9  are  illustrated  in  Fig.  2  over 
the  same  wide  range  of  </>,  po  and  Tq  as  in  Fig.  1.  Remarkably,  both  types  of  plots  display  a  quasi-linearity. 
For  H2O,  eventually,  an  asymptotic  behavior  is  seen  at  0  ~  0.6  indicative  of  the  reaction  having  reached 
completion,  consistent  with  the  information  from  Fig.  1. 

The  above  findings  suggest  that  instead  of  computing  Nc  from  a  rate  equation  (i.e.  as  a  process  variable), 
one  may  simply  fit  the  Nc/ {<j>xN*)  curves  shown  in  Fig.l;  this  fitting  must  be  performed  in  the  four  parameter 
space  {To,(j),N* ,T).  Similarly,  instead  of  considering  O2  and  H2O  as  progress  variables,  it  seems  reasonable 
to  functionally  fit  the  slope  of  the  curves  in  Fig.  2  and  use  this  fit  in  routine  calculations  to  predict  these 
two  major  species.  Thus,  this  examination  of  our  conceptual  model  for  n-heptane  using  the  LLNL  skeletal 
kinetics  in  CHEMKIN  II  leads  to  further  reducing  the  number  of  species-related  process  variables  from  12 
to  9.  These  9  process  variables  are  the  molar  densities  of  the  lights:  CO2,  CO,  H,  H2,  CH4,  II2O2,  C2II2, 
C2H4,  CH2O. 


B.  Assessment  of  the  concept’s  predictive  capabilities  using  an  ideal  model 

Before  launching  into  the  task  of  performing  functional  fits  for  rates  and  Xi,  we  deemed  it  instructive  to 
assess  the  predictive  capability  of  the  concept  by  utilizing  ideal  functional  fits  extracted  from  the  LLNL 
skeletal  mechanism  using  CHEMKIN  11.  These  ideal  functions  were  obtained  in  tabular  form,  every  5°K, 
and  served  as  input  to  our  conceptual  model.  Examples  of  the  results  are  portrayed  in  Fig.  3  for  both  Nc 
and  T  profiles;  indeed,  predicting  the  energetics  and  more  particularly  the  ignition  time,  tign,  is  an  important 
goal  of  the  model. 

Unsurprizingly  according  to  Fig.  1,  Fig.  3a  shows  that  Nc  is  excellently  predicted  by  our  conceptual 
model.  However,  whereas  the  value  of  the  combustion  T  is  also  excellently  predicted,  the  predicted  ignition 
time  is  longer  than  that  of  the  skeletal  mechanism.  This  epitomizes  the  missing  information  resulting  from 
the  reduction  of  the  LLNL  skeletal  mechanism  from  160  to  9  progress  variables.  Physically,  the  missing 
information  is  that  of  how  the  heavy  species  evolve  over  the  relatively  long  initial  time  of  nearly  constant 
T.  During  this  time,  although  T  is  almost  constant,  heavy  species  decomposition  provides  the  radicals 
responsible  for  eventual  ignition.  Experimentally,  it  is  observed^®  that  is  an  extremely  sensitive  quantity 
function  of  Tq  and  even  a  degree  K  change  in  To  makes  approximately  10%  change  in  tign-  Thus,  the  accurate 
prediction  of  is  a  very  stringent  test  for  a  model.  From  the  modeling  viewpoint,  this  means  that  the 
neglected  scales  of  the  initial  pre-ignition  process  must  be  reintroduced  through  very  careful  modeling.  This 
modeling  is  described  next. 

C.  Functional  fits  to  emulate  the  LLNL  skeletal  mechanism 

The  functional  fits  necessary  to  complete  the  model  fall  into  two  categories:  fits  for  the  rates  Kh,KGi, 
RLi  =  KLi/KGi,  for  N^  and  for  the  Xi  of  the  quasi-steady  light  species;  and  fits  for  the  slopes  of  the 
unsteady  quantities  No2,  and  Nh2o- 

1.  Fits  for  the  rate  quantities  and  quasi-steady  light  speeies  molar  fractions 

According  to  Eq.  8,  to  recover  the  value  of  T  in  the  reduction  scheme,  the  heavy  species  model  should  focus 
on  an  appropriate  representation  of  Gp^u  and  Kh-  Plots  (not  shown)  of  Gp^t  at  various  po  values,  calculated 
using  the  LLNL  model  over  a  wide  range  of  </>,  show  that  these  curves  exhibit  a  very  modest  variation  over 
the  range  of  strong  Nc  decay;  moderate  (i.e.  up  to  50%)  variation  occurs  only  after  Nc  values  are  very 
small  to  negligible.  This  indicates  that  the  T  recovery  is  primarily  governed  by  Kh  as  far  as  heavy  species 
modeling  is  concerned. 

The  model  for  Kh,KGi,  RLi  and  Xi  is  constructed  considering  these  quantities  as  functions  of  T  or  9, 
with  N*,4>  and  Tq  being  parameters.  Quantities  Kh,  KGi  and  Xi  are  split  into  a  very  low  T  region,  i.e. 
9  ^  10^^,  termed  the  incubation  region  in  which  there  are  very  slow  changes;  a  low- rate  portion  corresponding 
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to  10^^  ^  0  ^  0.2  that  exhibits  a  maximum,  K^x,  and  a  minimum,  Kmn]  and  a  high  T,  9  ^  0.2,  high-rate 
portion  termed  the  fast-rate  region  in  which  a  continuous  increase  is  commonly  observed.  Conversely,  RLi 
first  exhibits  a  minimum,  and  then  a  maximum.  Thus,  the  functional  fits  are  performed  in  three  separate 
regions,  with  connection  constraints  between  consecutive  regions  (the  ultimate  fits  are  piecewise  continuous). 
Details  of  the  fitting  procedure  are  given  in  Harstad  and  Bellan.^^ 

The  need  to  introduce  a  model  to  capture  tign,  as  explained  in  section  B,  as  well  as  an  extreme  sensitivity 
of  Kh  on  To  (e.g.  Kh  ~  Tg^  for  9  ~  10^^)  meant  that  special  care  should  be  devoted  to  reproduce 
this  dependency.  Thus,  we  developed  functional  fits  to  the  slopes  dKh/dT  at  Tg  and  made  corresponding 
adjustments  to  the  rates  at  0  <  10^^.  These  turned  out  to  be  insufficient  to  reproduce  tigm  and  an  adjustment 
was  introduced  to  the  initial  slope  for  (for  9  <  10^^)  by  using  a  multiplying  factor  determined  by  trial- 
and-error  calculations,  so  as  to  best  match  the  reduced  model  predictions  compared  to  those  of  the  skeletal 
mechanism. 

This  model  was  used  for  Kh,  for  the  quasi-steady  gain  rates  KGi  (where  i  stands  for  CO2,  H,  CO,  H2, 
CH4,  H2O2,  C2H2,  C2H4,  CH2O),  for  the  quasi-steady  loss  rate  RLi  (where  i  stands  for  H  and  H2,  this  rate 
being  negligible  for  the  other  light  species)  and  for  the  quasi-steady  light  species  molar  fractions  Xi  (O,  CH, 
CH2  CHa,  HO,  HCO,  HO2,  HC2,  C2H3). 

Selected  plots  of  Kh  are  presented  in  Fig.  4;  KGh2  and  KGh  are  displayed  in  Fig.  5  and  the  ratio  RLi 
is  correspondingly  illustrated  in  Fig.  6  for  f  =  H  and  f  =  H2.  The  rates  Kh  exhibit  a  variation  by  as  much 
as  seven  orders  of  magnitude,  and  the  task  of  developing  curve  fits  over  this  extended  regime  in  the  four 
parameter  space  (Tg,  0,  N*,9)  is  far  from  trivial.  Despite  the  difficulty  of  developing  accurate  curve  fits,  the 
agreement  between  the  fits  and  the  LLNL  skeletal  mechanism  in  excellent  to  good  over  pg  G  [10, 40]  bar, 
Tg  G  [600,800]  K  and  ^  G  [0.5,4]  and  even  extends  to  values  as  small  as  ^  =  1/4  (not  shown).  The  same 
comments  hold  for  KG  and  RL.  Discrepancies  between  fits  and  the  LLNL  skeletal  mechanism  past  9  ~  0.6 
are  unimportant  for  the  reduction  concept  since  Nc/Ncfi  0(1)  past  that  station. 

For  the  quasi-steady  species,  we  chose  to  display  No  because  of  the  high  reactivity  of  the  O  radical  and 
Noh  because  the  OH  species  is  generally  considered  as  indicative  of  the  flame  location.  The  quantities  Ni 
were  computed  as  Ni  —  XiN^  where  N^  =  '^i^ughts  Results  are  presented  in  Fig.  7  and  show  that  the 
agreement  between  fits  and  the  skeletal  mechanism  is  very  good,  including  some  extreme  (p  values  in  the 
lean  regime. 

2.  Fits  for  No2  and  Nh2o 

Two  quantities  were  defined  that  constrain  the  evolution  of  N02  and  Nh2o  '■  the  initial  value,  NI2,  and  the 
asymptotic  value  Nf^2o-  h'hs  for  N02  and  Nh2o  were  constructed  as  follows: 


N02INI2  =  max(1.0  -  Ao2  X  9,  0.0),  (12) 

Nh2o/Nh2o  =  aam(1.0,  Ah2o  X  &)■  (13) 

Fits  to  the  slopes  A02  and  Ah2o  were  constructed  as 

Ao2  =  Co2{<I>,To)  X  X  (iV*)0-05  (14) 

Ah2o  =  Ch2o{P,To)  X  {N*y  (15) 


where  e{(j))  =  0.2  x  p  x  exp(— 1.4(/)  <  0.0525.  For  computing  functions  C'o2(</>,  Tg)  and  C'?i2o(</>,  Tg)  a  set  of  f 
values  is  first  chosen.  For  a  specified  value  in  the  set,  Go2{4',To)  and  Gh2o{4',TQ)  are  fitted  as  a  polynomial 
function  (up  to  quadratic)  of  To/Tref-  For  arbitrary  p,  an  interpolation  is  performed  over  p  using  C02  and 
Ch2o  values  at  the  specified  f  values  in  the  chosen  set. 

Figure  8  shows  an  example  of  the  obtained  fits  versus  the  LLNL  skeletal  mechanism.  The  agreement 
ranges  from  excellent  to  good  and  is  typical  of  that  obtained  over  a  wide  range  of  (Tg,  (/>,  A^*). 

D.  Model  predictions 

The  model  predictions  consist  of  the  timewise  temperature  evolution  T{t),  Ugn  as  extracted  from  the  T{t) 
profile,  and  the  timewise  evolution  of  the  unsteady  light  species  with  the  exception  of  O2  and  H2O  which 
are  modeled  as  shown  in  section  2. 
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1.  Temperature  profiles 

Displayed  in  Fig.  9  are  the  reduced  model  predictions  compared  to  those  of  the  skeletal  mechanism,  spanning 
a  wide  range  of  </>  values  and  two  values  of  each  Tq  and  Pq.  Although  tign  is  well  predicted  for  all  cases, 
the  ultimate  combustion  temperature  is  not  well  predicted  at  very  lean  conditions  (</>  =  1/4)  or  relatively 
low  pressure  (i.e.  po  =  10  bar).  It  seems  likely  that  with  additional  modeling  effort  these  limitations  could 
be  overcome,  but  we  consider  that  at  this  ‘proof  of  principle’  stage,  the  additional  effort  is  not  warranted. 
On  the  favorable  side,  at  the  high  temperature  conditions  representative  of  diesel,  gas  turbine  and  HCCI 
engines,  the  model  performs  very  well;  the  same  holds  for  the  very  rich  conditions  (e.g.  </>  =  4)  at  which  soot 
forms. 

2.  Ignition  times 

The  very  successful  use  of  the  adjusted  initial  Kh  slopes  instrumental  in  predicting  tign  is  shown  in  Fig.  10. 
Noteworthy,  the  ordinate  axis  is  the  typical  logarithmic  one^*  for  Figs.  10a  and  10b,  but  it  is  here  linear 
for  Fig.  10c,  so  deviations  between  reduced  model  and  skeletal  mechanism  observed  in  Fig.  10c  are  indeed 
logarithmically  very  small.  The  predictions  in  Fig.  10  are  excellent  over  a  wide  range  of  Tq  including  as  low 
a  temperature  as  600  K,  over  a  range  of  </>  including  mixtures  as  rich  as  </>  =  4,  and  as  lean  as  4>  =  1/3  and 
Pq  from  10  to  40  bar. 

3.  Predicted  unsteady  species 

The  prediction  of  several  unsteady  light  species  is  illustrated  in  Figs.  11  and  12.  These  are  results  obtained 
from  solving  Eq.  6  with  the  modeled  TZi.  The  Ni{0)  is  computed  from  using  Ni{t)  and  T{t),  which  allows 
computing  9{t). 

We  chose  to  display  and  Nh2  because  they  both  depend  on  the  KG  and  RL  fits  and  thus  they  provide 
a  stringent  test  of  the  model  given  the  possibility  of  additive  error.  A/o  is  chosen  since  CO  is  one  important 
product  of  incomplete  combustion.  N^hA  is  selected  as  a  representative  intermediary  in  the  oxidation  process. 
As  discussed  above,  only  prediction  up  to  0  ~  0.6  should  be  considered,  as  past  that  station  the  reaction  is 
practically  finished. 

Generally,  the  agreement  between  the  reduced  model  and  the  skeletal  mechanism  is  very  good,  including 
at  po  =  10  bar  for  which  the  combustion  temperature  was  not  well  predicted.  With  the  exception  of  Nh  for 
which  the  comparison  is  just  as  favorable  at  ^  =  4  as  at  </>  =  1/3,  for  all  other  species  the  lean  region  is  much 
better  reproduced  by  the  model  than  the  rich  region.  The  fact  that  good  agreement  is  obtained  for  such 
extreme  f  values  is  desirable  since,  as  already  stated  such  values  are  generally  not  considered  in  advanced 
reduced  schemes.  Also,  the  results  seem  to  be  slightly  more  accurate  with  increasing  Tq. 

IV.  Preliminary  assessment  of  the  n-heptane  model  for  iso-octane 

To  test  the  capability  of  the  proposed  technique  for  other  alkanes,  iso-octane  was  chosen  as  the  next  one 
in  degree  of  difficulty.  Not  only  is  iso-octane  the  next  higher  C  number  species  with  respect  to  n-heptane, 
by  also  these  two  species  are  prominently  featured  in  JP-8  surrogate  fuels  because  they  both  display  the 
NTC  behavior  of  the  heavier  real  fuels  in  the  ignition  regime  and  the  feature  of  two-stage  ignition. The 
level  of  modeling  difficulty  increases  substantially  for  the  n-heptane  mechanism  having  160  species  (progress 
variables)  and  1540  reactions  to  the  iso-octane  full  (rather  than  skeletal)  oxidation  mechanism^®  containing 
857  species  and  3606  reactions.  The  substantial  effort  of  decomposing  in  constituents  the  heavy  species  of 
the  857-species  ensemble  led  to  the  iso-octane  ensemble  of  constituents  being  essentially  the  same  as  for 
n-heptane,  with  the  sole  addition  of  the  carbon  atom  C.  That  is,  now 

14 

N,  =  J2Nk  (16) 

k=l 

and  the  ensemble  of  constituents  is  (CHa,  CH3,  CH,  C2H3,  C2H2,  C2,  HC2,  CO  (keto),  HCO,  HO,  HO2, 
00,  O,  C). 

There  are  four  questions  that  are  of  interest  is  assessing  the  n-heptane  conceptual  model  for  iso-octane: 
(1)  Is  0  defined  in  Eq.  3  still  a  valid  similarity  variable  in  that  Nc/{(j)  x  N*)  as  a  function  of  9,  such 
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as  shown  in  Fig.  1  for  n-heptane?  (2)  If  the  similarity  variable  holds,  are  the  concentrations  of  O2  and 
H2O  decreasing  and  increasing,  respectively,  in  a  quasi- linear  manner  with  0?  (3)  Is  the  ensemble  of  light 
species  partitionable  into  the  same  two  sub-ensembles  of  quasi-steady  and  unsteady  species,  and  do  these 
sub-ensembles  contain  the  same  species  as  for  n-heptane?  (4)  If  the  three  above  items  hold,  are  the  functional 
fits  derived  for  n-heptane  also  valid  for  iso-octane?  The  preliminary  results  derived  for  iso-octane  allow  to 
cautiously  affirmative  answer  the  first  two  questions,  with  the  understanding  that  a  much  more  extensive 
study  is  necessary  for  incontrovertible  affirmation. 

Figure  13a  illustrates  Nc/{4'  ^  ^*)  as  a  function  of  9  at  three  ^  values  and  two  Tq  values.  We  first 
note  that  the  temperature  normalization  holds  in  that  Nc/{<i)  x  N*)<^  1  by  0  ~  0.6.  Second,  all  plots  of 
Nc/{4'  X  ^*)  approximately  coincide,  showing  for  this  limited  set  of  cases  that  the  similarity  concept  holds. 
Displayed  in  Fig.  13b  are  the  temperature  profiles  corresponding  to  the  Fig.  13a  cases,  showing  that  the 
choice  of  the  runs  spans  more  than  one  order  of  magnitude  range  in  tign- 

The  O2  and  H2O  concentrations  versus  6  are  plotted  in  Figs.  14a  and  14b.  The  quasi-linear  aspect  of  the 
curves  is  very  similar  to  that  observed  for  n-heptane  in  Fig.  2,  indicating  that  the  same  behavior  prevails. 

Although  considerably  more  tests  are  necessary  to  ascertain  that  the  concept  and  model  developed  for 
n-heptane  also  holds  for  iso-octane,  the  preliminary  results  seem  promising. 

V.  Conclusions  and  discussion 

A  kinetic  reduction  model  has  been  developed  based  on  the  concept  of  constituents  representing  the 
evolution  of  all  heavy  species,  and  on  light  species  representing  the  complementary  chemistry  to  that  of 
heavy  species.  The  light  species  fall  into  two  categories:  quasi-steady,  for  which  no  differential  equation 
must  be  solved,  and  unsteady  light  species  some  of  which  are  progress  variables  and  others  which  have 
fitted  rates  of  their  evolution  as  function  of  the  similarity  variable.  A  thorough  examination  of  the  LLNL 
skeletal  mechanism  for  n-heptane  revealed  that  it  is  possible  to  empirically  define  a  similarity  variable  which  is 
function  of  the  initial  temperature,  initial  pressure  and  equivalence  ratio,  and  for  which  a  non-dimensionalized 
global-constituents’  molar  density  exhibits  a  self-similar  behavior  across  initial  temperatures,  initial  pressures 
and  equivalence  ratios.  The  significance  of  this  finding  is  that  there  is  a  reduction  in  the  dimensionality  of 
the  problem.  We  found  a  lower-dimensional  manifold  which  is  not  to  be  confused  with  ILDM  because  the 
latter  is  for  species,  whereas  the  present  model  is  for  W-  The  fact  that  the  constituents  rather  than  the 
heavy  species  are  important  can  be  directly  traced  to  the  fact  that  the  heavy  species  decompose  and  it  is 
the  reaction  of  these  products  of  decomposition  that  matters  for  the  energetics.  Thus,  rather  than  solving 
a  differential  equation  to  obtain  the  evolution  of  the  global  constituents,  their  behavior  may  be  functionally 
fitted,  thereby  reducing  the  number  of  progress  variables.  Further  examination  of  the  skeletal  mechanism 
revealed  that  the  molar  density  of  oxygen  and  that  of  water  displayed  a  quasi-linear  variation  with  the 
similarity  variable  over  the  entire  range  of  parameters  surveyed.  The  suggestion  was  that  the  molar  densities 
of  these  two  unsteady  light  species  could  be  functionally  modeled  rather  than  computed  from  differential 
equations. 

The  concept  was  tested  by  using  the  n-heptane  LLNL  skeletal  mechanism  in  conjunction  with  CHEMKIN 
11.  That  is,  instead  of  a  model  consisting  of  functional  fits,  we  used  an  ideal  model  represented  by  tables 
extracted  from  the  LLNL  skeletal  mechanism  in  the  form  needed  in  our  conceptual  model.  Our  conceptual 
model  was  found  to  be  sound,  but  the  ignition  time  it  predicted  was  too  long  compared  to  the  template. 
The  meaning  of  this  discrepancy  is  that  the  neglect  of  the  chemical  processes  during  the  initial  time  of  very 
small  temperature  changes  results  in  the  overprediction  of  the  ignition  time.  The  situation  is  equivalent 
to  turbulent  modeling  using  the  Large  Eddy  Simulation  (LES)  concept  in  which  the  spatial  filtering  of 
scales  results  in  the  need  to  re-introduce  the  effect  of  these  small  scales  through  models.  Following  the  LES 
modeling  philosophy,  a  small-scale  model  was  developed  and  implemented  at  temperatures  very  close  (1-2 
degrees  K)  to  the  initial  temperature,  to  be  used  with  the  conceptual  model. 

Moreover,  also  paralleling  the  LES  concept,  fits  in  functional  form  were  developed  for  the  sensible  enthalpy 
production  due  the  constituents’  reaction,  for  the  quasi-steady  light  species  mole  fractions  and  for  the  reaction 
rate  of  the  unsteady  light  species.  The  reaction  rate  of  each  unsteady  light  species  was  decomposed  into 
contributions  from  the  lights  which  were  directly  taken  from  the  LLNL  skeletal  mechanism  and  contributions 
from  the  heavy  species  which  were  modeled.  The  heavy  species  rate  contribution  was  further  categorized 
into  a  gain  rate  and  possibly  a  loss  rate,  which  were  individually  modeled.  The  modeling  effort  was  by  no 
means  trivial  as  the  rates  and  molar  fractions  varied  by  as  much  as  seven  orders  of  magnitude  and  this 
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variation  needed  to  be  captured  in  the  four  parameter  space  of  the  initial  pressure,  the  initial  temperature, 
the  equivalence  ratio  and  the  temperature. 

The  consequence  of  providing  functional  fits  over  a  wide  range  of  equivalence  ratios  encompassing  mixtures 
as  rich  as  ^  =  4  and  as  lean  as  ^  =  1/4  was  that  generally  the  fits  were  very  good  to  excellent,  but  sparse 
regions  of  good  to  fair  results  could  also  be  identified.  The  ignition  times  were  excellently  to  very  well 
reproduced  by  the  model,  and  so  were  some  major  species  (e.g.  O2  and  H2O).  Other  major  species  (e.g. 
CO)  and  OH,  which  is  an  indicator  of  the  flame  location,  were  also  well  reproduced,  but  the  value  of  the  hnal 
combustion  temperature  was  not  always  well  predicted.  To  improve  the  combustion  temperature  prediction, 
one  could  either  refine  the  functional  fit  or  use  the  ideal  model  provided  by  the  tables  extracted  from  the 
LLNL  mechanism  using  CHEMKIN  II  as  done  in  section  B,  but  now  in  conjunction  with  our  model  in  the 
vicinity  of  Tq  so  as  to  capture  the  ignition  time.  To  improve  light  species  predictions,  one  could  either  refine 
their  functional  fit  if  the  species  is  quasi-steady,  or  refine  the  functional  hts  for  KG  and  RL  if  the  species 
is  unsteady;  alternately,  one  could  use  the  the  ideal  model  provided  by  the  tables  extracted  from  the  LLNL 
mechanism  using  CHEMKIN  H. 

Instead  of  improving  the  n-heptane  model,  we  turned  attention  to  investigating  whether  the  proposed 
concept  also  holds  for  higher  C  alkanes,  such  as  iso-octane.  Preliminary  results  indicate  that  the  conceptual 
model  of  the  similarity  variable  and  the  quasi-linear  dependence  of  the  O2  and  H2O  concentrations  of  the 
similarity  variable,  indeed  holds;  although  it  remains  to  be  seen  if  the  present  curve  hts  are  valid  for  higher 
C  alkanes  or  must  be  modihed. 

Finally,  the  usefulness  of  our  model,  which  contains  only  nine  species-related  progress  variables,  is  that 
it  is  less  computational  intensive  than  other  reduced  models  having  0(50)  species-related  progress  variables 
(e.g.  Lu  and  Law^®)  while  giving  much  improved  predictions  compared  to  the  expedient  one-step  reaction 
models.  As  usual,  information  loss  increases  with  increasing  chemical  kinetics  simplihcation;  the  balance 
between  the  desired  information  to  be  obtained  from  exercising  the  model  and  computational  cost  is  the 
paramount  consideration.  Since  results  for  most  reduction  schemes  are  presented  in  a  much  more  restricted 
(j)  region  than  the  (j>  G  [1/4,4]  considered  here,  it  is  difficult  to  directly  evaluate  our  model  with  respect  to 
other  models.  It  is  understood  that  for  any  given  model,  predictions  will  improve  if  a  less  ambitious  validity 
regime  is  considered. 
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Mo/Ra  /i°  h‘^ 


Ha 

O2 

Na 

C 

HaO 

COa 

N 

H 

HO 

HOO 

O 

CO 

HCO 

CH4 

CH3 

CHa 

CH 

CaHa 

HCa 

Ca 

NO 

NOa 

HaOa 

HCOH 

CaHa 

CaH4 


0.0 

0.0 

0.0 

717 

-241.5 

-393.5 

473 

218.0 

38 

10.5 
249.2 
-110.5 

43.1 
-74.6 
146 
390 
596 
300 
566 
838* 
90.3 

33.2 
-136 
-109 
228 

52.5 


241.5 
0.0 
0.0 
nil 
0.0 
0.0 
473 
339 
159 
131 

249.2 
283 
558 
802 
902 
1025 
1110 
1449 
1474 
1625 
90.3 

33.2 
106 

526.5 
1257 
1323 


3.282 

3.476 

3.388 

2.50 

3.688 

4.690 
2.50 
2.50 
3.385 
4.150 
2.536 
3.426 
4.154 
3.797 
4.440 
3.973 
3.220 
5.1 
4.434 
4.58 
3.533 

4.691 
5.269 
4.27 
5.368 
5.383 


0.400 

0.5663 

0.469 

0.0 

1.217 

1.390 

0.0 

0.0 

0.3637 

1.307 

0.0 

0.4749 

1.2875 

4.305 

2.249 

1.3015 

0.7136 

3.5 

1.404 

0.0 

0.4508 

1.151 

1.880 

2.546 

2.294 

4.676 


*Alternate  value  of  832  from  the  CRC  tables. 


Table  1.  Thermodynamic  properties  of  molecules  and  free  radicals,  and  hP  (heats  of  formation  and  combus¬ 
tion,  respectively),  in  kJ/mol;  constants  a'*  and  for  molar  heat  capacity  in  the  form  Cp/Ru  =  -\-b^ln{T/Tref)'’, 

=  298.15  K.  "Mo"  denotes  "molecule".  "Ra"  denotes  "radical". 
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Figure  1.  Similarity  plots  of  parameter  Nc/{(p  x  N*)  versus  6  at  (a)  po  =  20  bar  and  (b)  0  =  1  using  the  LLNL 
skeletal  mechanism  in  conjunction  with  CHEMKIN  II.  Tq  is  in  K. 
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Figure  2.  Oxygen  and  water  molar  densities  versus  6  as  extracted  from  the  LLNL  skeletal  mechanism  in 
conjunction  with  CHEMKIN  II. 


Figure  3.  Molar  density  of  constituents  (a)  and  temperature  (b)  obtained  with  our  conceptual  model  and 
tabular  information  from  the  LLNL  skeletal  mechanism  using  CHEMKIN  II.  The  symbols  are  from  the  LLNL 
skeletal  mechanism  in  conjunction  with  CHEMKIN  !!(□□□  po  =10  atm;  O  O  O  P0=  40  atm),  and  the 
lines  are  obtained  using  tabular  information  from  the  LLNL  skeletal  mechanism  with  CHEMKIN  II  in  our 
conceptual  model  ( - po  =10  atm; - P0=  40  atm). 
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Figure  7.  Modeling  results  for  the  quasi-steady  molar  densities  No  and  Noh  and  comparison  with  those  of  the 
skeletal  mechanism  at  various  To,po  and  (p  conditions. 


16  of  21 


American  Institute  of  Aeronautics  and  Astronautics 


Figure  8.  No2  and  Nf^2o  versus  0  at  fixed  po  =  10  and  40  bar  for  <p  =  1  and  Tq  =  700  and  800  K. 


0.0025 


0.005 

t(s) 


0.0075 


0.01 


Figure  9.  T{t)  as  predicted  by  the  reduced  model  compared  to  the  LLNL  skeletal  mechanism  predictions  for 
a  variety  of  conditions. 
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Po  (bar) 


Figure  10.  tign  as  function  of  To,0  and  pQ.  Where  not  indicated  on  the  plots,  Tq  is  in  degrees  K  and  po  is  in 
bar. 
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Figure  11.  Predictions  of  the  reduced  model  for  and  Nh2  compared  to  those  of  the  skeletal  mechanism  at 
various  To,po  and  (p  conditions. 
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Figure  12.  Predictions  of  the  reduced  model  for  Nco  and  Nch4  compared  to  those  of  the  skeletal  mechanism  at 
various  To,po  and  (p  conditions. 
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Figure  13.  Iso-octane  results  obtained  using  the  LLNL  full  mechanism  in  conjunction  with  CHEMKIN  II:  (a) 
similarity  plot  of  parameter  Nc/{(p  x  A^*)  versus  0  at  po  =  40  bar  and  (b)  T  time-wise  evolution.  Tq  is  in  K. 


Figure  14.  Oxygen  and  water  molar  densities  versus  9  as  extracted  from  the  iso-octane  LLNL  full  mechanism 
in  conjunction  with  CHEMKIN  II.  po  is  in  bar,  Tq  is  in  K. 
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APPENDIX  4 


Computation  of  Laminar  Premixed  Flames  Using 
Reduced  Kinetics  Based  on  Constituents  and  Species 

Kenneth  G.  Harstad^’*and  Josette  Bellan^’**d 
Jet  Propulsion  Laboratory* ,  California  Institute  of  Technology, 

Pasadena  CA  91109-8099 

California  Institute  of  Technology** ,  Mechanical  Engineering,  Pasadena,  CA  91125 


A  model  is  proposed  for  quasi-one-dimensional  steady  flame  development  in  the  con¬ 
figuration  of  an  inviscid,  premixed  fuel  jet  injected  into  air.  The  governing  equations 
are  written  within  the  framework  of  a  reduced  kinetic  model  based  on  constituents  and 
species.  The  reduced  kinetic  model,  previously  exercised  in  a  constant-volume  perfectly- 
stirred  reactor  mode,  has  been  successful  at  predicting  ignition  and  combustion  product 
and  temperature  evolution  for  n-heptane,  iso-octane,  PRF  fuel  combinations,  and  mixtures 
of  iso-octane  with  either  n-pentane  or  iso-hexane.  The  differential  governing  equations  have 
the  option  of  an  axially  variable  area  and  they  are  coupled  with  a  real  gas  equation  of  state. 
The  flame  development  model  accounts  for  a  full  diffusion  matrix,  and  thermal  conductivity 
computed  for  the  species  mixture.  Preliminary  results  are  presented. 


I.  Introduction 

To  predict  the  initiation  and  evolution  of  turbulent  reacting  flows  in  practical  devices  one  must  couple 
chemical  kinetics  and  turbulence  models.  Because  codes  resulting  from  these  models  are  very  computationally 
intensive,  the  chemical  kinetic  mechanisms  cannot  be  utilized  in  their  entire  complexity  and  instead  they  are 
reduced  to  mechanisms  that  are  sought  to  be  compact  as  well  as  reliable.  Recently  a  reduced  kinetic  model 
relying  on  constituents  and  species^’ ^  has  been  developed,  which  was  successful  in  reproducing  ignition  time, 
and  combustion  product  and  temperature  evolution  for  n-heptane,  iso-octane,  PRF  fuel  combinations,  and 
mixtures  of  iso-octane  with  either  n-pentane  or  iso-hexane. 

In  this  conceptual  model,^’^  the  ensemble  of  species  in  a  full  (or  skeletal)  mechanism  is  partitioned  into 
heavies  (carbon  number,  n  ^  3)  and  lights  (the  remaining  of  the  set).  The  heavies  can  be  either  radicals  or 
stable  species.  The  lights  are  oxygen,  nitrogen,  the  hnal  combustion  products  and  light  radicals/molecules 
(e.g.  CH3,  CH4,  H2O2).  The  heavies  are  not  treated  as  a  species  set,  but  as  a  set  of  base  constituent  radicals 
which  correspond  to  radicals  as  used  in  group  additivity  theory.^  Sometimes  light  species  are  radicals  which 
may  have  the  same  chemical  formula  as  constituents,  however,  the  difference  between  constituents  and  these 
light  species  is  that  the  later  are  unbound  to  other  chemical  entities  whereas  the  former  are  bound  to  other 
chemical  entities,  i.e.  other  constituents.  For  each  constituent  k,  its  molar  density,  Nk,  is  the  sum,  over  all 
heavy  species,  of  the  count  of  the  constituent  in  each  heavy  species  multiplied  by  the  molar  density  of  that 
species.  The  constituents  are  not  just  based  on  atom  counts  and  their  element  composition  is  not  linearly 
independent,  however,  the  constituents  are  linearly  independent.  Thus,  the  constituents  are  independent 
structural  elements  which  have  individual  valence  bond  topologies. 

For  each  heavy  species,  and  thus  for  the  entire  set  of  heavy  species,  there  is  a  unique  and  complete  set  of 
constituents.  The  subset  of  constituents  that  have  a  very  small  contribution  to  the  total  constituent  count 
having  molar  density  N,,,  is  replaced  by  the  complementary  subset  of  constituents  already  accounted  for  in 
Nc-  The  rule  for  replacement  is  that  each  constituent  having  a  small  contribution  to  Nc  is  replaced  by  one 
having,  first,  close  valence  structure  and  second,  close  composition  to  the  constituent  neglected.  After  this 

*Senior  Engineer. 

^Senior  Research  Scientist,  AIAA  Fellow  (corresponding  author,  josette.bellan@jpl.nasa.gov). 
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replacement,  one  obtains  an  optimal  Nk  set.  It  is  this  optimal  set  of  constituents  that  is  used  in  the  kinetic 
reduction. 

To  test  the  capability  of  the  reduced  kinetic  model  in  the  a  simple  flow  model,  the  first  step  is  to  evaluate 
it  in  a  steady,  quasi-one-dimensional  configuration  before  undertaking  the  more  arduous  task  of  the  counter¬ 
flow  configuration  for  which  the  overwhelming  experimental  data  has  been  obtained.  This  first  step  is  the 
purpose  of  this  study.  In  this  paper,  we  first  present  the  governing  equations,  as  they  are  more  general 
than  the  typical  ones  for  steady  one-dimensional  flow  which  use  approximate  transport  properties  and  the 
perfect  gas  equation  of  state.  Further,  we  discuss  some  preliminary  results  obtained  with  the  model.  Finally, 
we  elaborate  on  the  future  work  necessary  before  such  a  model  is  ready  for  incorporation  in  turbulent  flow 
models  and  codes. 


II.  Governing  equations 

The  configuration  of  interest  is  that  of  a  steady,  round  jet  in  a  system  of  coordinates  {x,  r,  0)  where  it 
is  assumed  that  there  is  uniformity  in  the  0  direction  and  thus  the  coordinates  of  interest  are  {x,r).  The 
entrance  of  the  jet  {x  =  0)  is  labeled  by  the  subscript  0  and  the  far  field  boundary  is  at  a:  =  L.  Within  the 
quasi-one-dimensional  framework,  all  quantities  are  averaged  in  the  r  direction,  and  thus  the  jet  boundary 
is  located  at  r  =  R{x).  Therefore,  the  jet  area  is  A{x)  =  'xR^{x).  A  quantity  ^{x)  =  Ro/R{x)  is  defined  to 
track  the  jet  spacewise  variation. 


A.  Continuity  and  momentum 


We  define  the  mass  density  p,  and  the  axial  and  radial  velocities,  u  and  v.  Radial  averages  are  designated 
by  0,  as  for  example 


p{x) 


R 

^  j  P{x,r)  r  dr. 
0 


(1) 


In  the  entire  derivation  below,  all  quantities  will  be  considered  as  radial  averages,  but  the  ()  symbol  will  be 
omitted  for  notation  simplicity. 


1.  Continuity  equation 

•  • 

Defining  the  mass  flow  rate  M,  the  steady  flow  assumption  implies  puA  =  M  =  constant,  which  means  that 


pu  =  M/A  =  M/wR^ix). 


Replacing  eq.  2  in  the  steady  continuity  equation 


1  djrpv)  d{pu)  ^ 

r  dr  dx 

leads  to 

M  dR 


(2) 

(3) 

(4) 


2.  Momentum  equation 

Under  the  assumption  of  an  inviscid  fluid,  the  steady  momentum  equation  components  are 


f  d  d\ 

/a  a  \ 


dp 

dx 

dp 


0 

0. 


(5) 

(6) 


To  make  progress  in  obtaining  an  equation  for  the  jet  area,  we  expand  p  and  p  ^  in  even  powers  of  r  (due  to 
symmetry  conditions),  as  a  means  of  taking  advantage  of  the  configuration  symmetry.  Thus,  we  postulate 
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that 


^  '^max  ^2n 

P  =  P0  +  ^Poul 

n—0 


~  7?2" 


r.2n 


n=0 


i?2ri 


where  cr„  and  /„  are  arbitrary  functions.  Replacing  these  expressions  in  eqs.  5  and  6  yields 


1  dan  _  . 

2  dx 


2diefn) 

dx 


.  d^  d 


d( 


dx 


dx 


n—1 


dx 


for  n  >  0. 


Combining  eqs.  9  and  10  for  n  =  1,  leads  to 

.4rf(eVl) 


dx 


dx 


f£|A^ 


dx 


dx 


=  0. 


(7) 

(8) 

(9) 

(10) 


(11) 


An  analysis  of  eq.  11  provides  insights  from  which  further  simplifications  are  obtained.  For  example,  near  a 
stagnation  point,  which  is  obtained  for  5^0,  the  second  term  on  the  left  hand  side  of  eq.  11  dominates  the 
first  term  which  can  then  be  neglected.  Also,  near  the  a:  =  0  boundary,  the  jet  is  nearly  radially  uniform, 
which  means  that  fi  =  0.  These  simple  considerations  indicate  that  /i  can  be  neglected  over  the  entire 
length  of  the  jet.  More  generally,  if  a;  7^  Xa,  where  Xa  is  one  of  the  elements  in  the  set  defined  by  those  x 


values  for  which 


^  1,  the  jet  area  can  be  computed  from 


e  =  l-(l-^J 


(12) 


where  is  the  eigenvalue  corresponding  to  the  variable  A(a:)  situation.  For  A  constant,  none  of  these  con¬ 
siderations  apply  and  the  momentum  equation  becomes  a  vehicle  to  compute  Vp,  but  under  the  assumption 
of  very  small  Mach  number  the  pressure  is  assumed  constant  for  the  calculation  of  the  equation  of  state 
(EOS). 


B.  Species  and  energy  equations 

The  complete  formulation  for  the  species  and  the  energy  equations  has  been  presented  elsewhere. However, 
that  formulation  was  for  species  and  it  is  not  applicable  here  where  the  model  is  for  both  constituents  and 
species.  In  this  reduced  kinetic  model,  the  evolution  of  the  global  constituents  molar  density  is  found  from 
input  tables  and  it  is  only  the  light-species  equations  that  are  progress  variables,  together  with  the  energy; 
however,  the  constituents  make  contributions  to  the  light,  and  these  contributions  must  be  now  formulated 
in  the  species  equations  in  terms  of  the  constituents.  To  this  end,  we  assign  indices  i,j  and  q  to  an  ensemble 
L  denoting  the  lights,  k  to  the  ensemble  C  of  the  constituents  and  l,m  to  the  ensemble  H  of  the  heavies. 
Ensemble  L  is  further  partitioned  into  that  of  the  quasi-steady  species,  LQS,  and  that  of  the  unsteady 
species,  LUNS.  All  quantities  are  assumed  averaged  over  the  r  direction,  so  there  is  only  a  variation  with  x. 


1.  Species  equation 


We  partition  the  species  molar  flux  in  the  light-species  governing  equation  in  two  categories  as  Jijight  which 
represents  the  contributions  on  light  species  i  of  all  lights,  and  Ji^heavy  which  represents  contributions  on 
light  species  i  from  all  heavies 


Ji  —  Ji, light  Ji, heavy 

(13) 

/  __  _  d\nT  dXj\ 

\  ^  / 

(14) 

—  Tl  —  1 

heavy  —  ^  ^  ^il  H,i 

(15) 
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where  n  is  the  molar  density,  X  is  the  mole  fraction,  T  is  the  temperature  and  quantities  Vij  are  pairwise 
diffusion  coefficients  for  the  lights  that  are  computable  from  the  binary  diffusion  coefficients  through  mixing 
rules;  ways  to  compute  binary  diffusion  coefficients  have  been  given  by  Harstad  and  Bellan®  and  mixing  rules 
were  derived  by  Harstad  and  Bellan.®  The  term  proportional  to  VT  includes  the  entire  Soret  effect  (in  Ji 
the  term  proportional  to  Vp  is  considered  negligible  because  it  is  usually  small  and  we  additionally  make  the 
small  Mach  number  approximation)  and  it  is  included  in  Ji^Hght  because,  as  for  Vij,  Dx^i  can  be  computed 
using  mixing  rules.  These  mixing  rules  generally  include  contributions  from  the  heavies,  however,  since  the 
individual  heavies  are  not  individually  tracked,  we  use  a  single  mean  heavy  species  dynamically  computed 
from  the  injected  alkane  fuel  as  it  pyrolyses.  The  goal  is  now  to  express  JH,i  as  a  function  of  a  constituent 
equivalent  of  a  mole  fraction. 

If  one  considers  a  constituent  Nk,  then 

Nk  =  Y,  CuNi,  N,  =  Y,Nk  (16) 

l  k 

where  Cki  is  the  count  of  constituent  k  in  species  1.  To  find  a  way  of  writing  the  heavies  as  a  function  of  the 
constituents,  one  must  invert  the  matrix  of  elements  Cki-  However,  because  there  are  fewer  equations  (i.e. 
constituents)  than  unknowns  (i.e.  heavies),  the  inversion  can  only  be  approximately  performed.  To  perform 
this  inversion,  we  select  the  Householder  transformation  which  when  utilized,  yields 

Xfn  —  ^  ^  CmkXk  where  Sml  —  ^  ^  CmkCkl  (11^) 

k  k 

with  5mi  being  the  Kronecker  symbol.  One  may  interpret  Nm  as  the  minimum  norm  of  heavy  molar  densities 
for  given  Nk-  This  norm  favors  heaviest  species  and  it  is  thus  most  accurate  at  early  times,  when  diffusion 
is  important.  Since  Xi  —  Ni/n, 


J, 


HA 


cHn(n) 

dx 


(18) 


Since  ultimately  only  ^k  is  available  in  the  model,  an  additional  transformation  is  needed  in  the 

expression  of  eq.  18.  Under  the  verified  assumption  that  individual  dominant  constituents  of  molar  density 
Nk  are  quasi-steady,  it  was  postulated^  that  there  exists  a  mole  fraction  XCk  such  that  Nk  —  NcXCk- 
Therefore, 


Y^VaCik  XCk  = 

Lk  I 


(19) 


Finally,  it  is  necessary  to  compute,  from  the  Vu  coefficients,  a  global  diffusion  coefficient  for  the  term 
accounting  in  the  formulation  for  the  global  constituent  molar  density.  For  this  purpose,  we  define 


V 


_  1 
HA  =  — 

n 


^  AilV;  =  ''^ViiXi  =  XhT^iH 
I  I 


(20) 


where  Xh  is  the  global  mole  fraction  of  the  heavies  and  ViH  is  a  mean  diffusion  coefficient  from  the  light 
species  i  to  the  ensemble  of  the  heavies.  If  one  defines  a  mean  total  constituent  count  as 


N, 


Cave=  —  ^y^CklXu 
n 


then  one  may  approximate  JH,i  by 


Jha  =  —NDh, 


This  leads  to 


Ji  =  -n 


k,l 

dlnCg 

dx 

dlnCg. 
dx 


H- 


■X,Dt, 


dlnT 
'  dx 


(21) 

(22) 

(23) 
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To  write  Ji  as  a  function  of  mass  fraction  gradients  and  take  advantage  of  the  partition  of  the  lights  set  into 
unsteady  and  quasi-steady  species,  we  now  assign  indices  i  and  j  to  the  unsteady  light  species  and  q  to  the 
quasi-steady  light  species.  Thus  the  total  unsteady  mole  fraction  is 


X„  =  l-X 


and  the  mixture  mean  mass  becomes 


H 

q 

X„, 


m  — 


y  21 

/  V  irii 


(24) 


(25) 


Introducing  these  definitions  in  eq.  23  yields 

Ji  = 


~n 


=  —n 


JXa  X — ^  ^  dXf, 


dx 


-x„n 


dlnCg 
dx 


X,D 


dlnT 


Vf 


dXg 

dx 


X - ^  \  dX^  ^ - T  dXf 


dx 


T,i- 


x„n 


dx 
d\nC, 


dx 


-  XiOxi- 


dliiT 

dx 


(26) 


having  used  dXj  =  {m/mj)dYj  -I-  Xjd\n{m)  and  having  defined  a  mean  unsteady-lights  mass  diffusion 
coefficient 

(27) 

^  i 

Writing  eq.  26  consistent  with  the  kinetic  reduction  modek’^  in  which  variables  Xq^X^  and  Cave  are  ob¬ 
tained  from  the  reduced  model  as  functions  of  the  temperature,  so  that,  for  example  dXjf /da;  =  {6XH/ST){dT/dx), 
leads  to 

dr] 


Ji  =  -n 


X — ^  .-r^n\  JYj 


rrij  dx 


B, 


where 


^X,^  +  y  (A,  -  +  XhV,h- 


6Xa 


dx 
din  a 


5T 


ST 


ave  2^0 


*  6T 


(28) 


(29) 


acts  as  an  effective  Soret  coefficient  for  the  reduced  model. 

Denoting  the  total  reaction  rate  of  species  z  by  7?.^,  the  governing  equation  for  species  of  mass  fraction 
P  is 

/  \ 

dx 


4-miT^i  -  J  i  4-mP 

M  , 


(30) 


In  accordance  with  the  reduced  kinetic  model 

TZi  =  TZi 


lights 


TZi 


heavies 


(31) 


where  is  computed  using  CHEMKIN  II  in  conjunction  with  the  LLNL  rate  data,  and  = 

NcKnet,i  where  Knet,i  is  computed  by  interpolation  from  tables  where  it  is  listed  as  a  function  of  T  for  fixed 
{To,pq,P  with  (j)  being  the  equivalence  ratio  and  subscript  0  denoting  initial  conditions. 


2.  Energy  equation 

We  write  the  energy  equation  under  the  assumption  that  the  Dufour  effect  has  a  negligible  contribution  to 
the  heat  flux.  Additionally,  since  in  our  model  the  heavy  species  are  combined  into  a  single  constituent, 
we  approximate  the  molar  enthalpy  of  all  heavy  species  averaged  over  the  molar  fluxes  by  the  enthalpy 
hn  —  (1/kfff)  which  is  the  mean  heavy-species  enthalpy  based  on  an  average  of  heavy  species  molar 

enthalpies  hi.  The  result  is 

9  ^ (32) 
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where  A  is  the  mixture  thermal  conductivity.  The  energy  equation  is  thus 


CpdT 


1  Ax—  \  _  A  J  ^ 

m  dx  dx  dx  I 


/  Cp,i  Cj 
^  '  rtii  niH 


p,H 


{rriiJi)  +  '^2/  —  .Rii7r.e/A^c.fi"ff  j" 


(33) 


where  Cp  is  the  heat  capacity  at  constant  pressure  and  Kh  is  computed  by  interpolation  from  tables^  where 
it  is  listed  as  a  function  of  T  for  fixed  {To,po,(l))- 


C.  Equation  of  state 

The  pressure  is  calculated  from  the  Peng-Robinson  (PR)  EOS 

_  RuT  Ojfnix 

{vpR  —  h  mix')  i^^PR  T  ^mix) 

where  vpR  is  the  molar  volume  and  Umix  and  bmix  are  functions  of  T  and  Xi  (see  Appendix  A). 


(34) 


D.  Transport  properties 

To  derive  transport  properties  consistent  with  the  kinetic  model  approach,  we  define  an  ensemble  of  all 
light  species  and  one  mole- averaged  heavy  species:  LU  H.  Indices  p,  n  and  r  refer  to  any  species  from 
this  ensemble.  The  transport  properties  under  consideration  are  the  diffusion  coefficients  the  thermal 
diffusion  factors  Dpp,  and  the  thermal  conductivity;  the  viscosity  does  not  enter  the  calculations,  since  an 
inviscid  situation  is  assumed,  however,  because  transport  property  calculations  of  thermal  conductivity  and 
viscosity  are  very  much  related,  we  borrow  from  methods  to  compute  viscosity  in  order  to  compute  the 
thermal  conductivity. 

For  the  diffusion  coefficients,  the  first  task  is  to  compute  the  binary  diffusion  coefficients  which  are  the 
building  blocks  of  the  pairwise  diffusion  coefficients.  To  this  end,  we  adopt  the  method  of  Harstad  and 
Bellan®  which  gives  (in  cgs  units) 


nD 


pn 


=  2.81  X 


fniT) 

2/3 


(35) 


where  the  subscript  C  denotes  the  critical  state,  /d(T)  =  {TredY  with  Ins  =  X)a=o  where  the 

a®  vector  has  elements  {-0.84211,  -0.32643,  -0.10053,  0.07747,  0.0127,  -0.00995},  and  rp  is  a  constant  0(1) 
which  provides  an  empirical  adjustment  for  the  specifics  of  the  collisional  interactions  of  a  selected  pair  of 
species.  Values  of  td  are  listed  in  Harstad  and  Bellan®  for  species  pairs  relevant  to  combustion.  The  heavy 
species  molar  mass  is  part  of  the  reduced  model  tables  and  decreases  as  the  reaction  proceeds. 

From  the  Dpn  ’s,  the  pairwise  diffusion  coefficients®  can  now  be  computed  using,  as  an  approximation  (a 
truncated  series  proposed  by  Ern  and  Giovangigli^ ) , 


V  — 

^pr  — 


-Dpn  — 


Xp  ^  ^  DpnCXp  ,nr  1 
n 

(i  +  ’^p) 


T)*T)* 


N 


^pUpn 

yip 


pn)  „  (f^p^p  +  fT'rtR’n)  +  ^^(VrfT'rR’r)) 


VI  = 


/  -^pn 


N 


vt 


D  ’ 

,  -^pn 

n^p 


where  the  mass  diffusion  factors  are  computed  as 


_  dXp  ainyp 

OiD.pn  —  +  ^p 


dX„ 


(36) 

(37) 

(38) 

(39) 

(40) 
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from  the  EOS,  with  7^  =  where  ^  is  the  fugacity  coefficient  and  the  superscript  o  denotes  the  pure 

{Xp  =  1)  limit. 

According  to  Harstad  and  Bellan® 


n—1  \  r=l 
\r^n 

T  _ 


(mr  +  mn)Dr 


(41) 


(42) 


where  Xp  is  the  species  p  thermal  conductivity  computed  according  to  eq.  50  and  is  a  weighting  factor 
given  by  eq.  45.  Factors  Cr  ra(^red,r  n)  are  closely  related  to  the  temperature  derivative  of  Dp„  taken  at 
constant  pressure,®  and  they  are  computed  from 


Cr  n{'^red,r  n)  —  ^  g 


n  +  Tred,  r  n  ^n{Tred  ,r  71) 


dSr 


dT, 


red', 


where  Tred,r  n  is  defined  in  Appendix  A. 

To  compute  A,  we  use  the  Wassiljewa-Mason-Saxona  method  described  in  Reid  et  al. 

A=  ^  X„c^QA„ 

nGLU// 

where  the  weighting  factors  oj^  are  computed  as  recommended  in  Reid  et  ah,® 

(CCQ)-I  =  J2^PnXn, 

3 

/  /v  \  1/2-1  2 

‘  +  (s;)  (I:) 


\/8(l  +  Sf) 


(43) 


(44) 


(45) 


(46) 


where  rjp  represents  the  viscosity  of  species  p.  To  compute  the  individual  species  viscosities,  we  adopt  a 
method  explained  in  Reid  et  al.*  whereby 


Vp  Vref,p  f R—T^3^red^p^  1 

where  the  subscript  R  —  T  stands  for  “Roy-Thodos”  and 

fR^T{Tred,p)  =  2.25(exp(0.0464Tred,p)  -  exp(-0.2412r,.ed_p)), 

2/3 


Vref,p  =  1-08  X  10  ■^(mprc,p)^/^ 

\  P’C,p  J 


(47) 

(48) 

(49) 


The  other  ingredient  entering  the  A  calculation  is  the  expression  for  Xp  which  is  here  computed  (in  cgs  units) 
using  the  Stiel-Thodos  method  (Reid  et  al.®)  as 


Xp  —  Xj^^f  p  \f R—T^red^p)  ^  ^C,p  ' f EjpiPred^p)^  , 


fx  =  3.75 


Cp  5 

Rp  2 


0.7862  -  0.71090  +  1.316802 


(50) 


(51) 


where  fx  is  a  function  of  the  acentric  factor  0  given  by  the  Chung  at  al.  formula,®  fR~T{Tred,p)  of  eq.  50 
accounts  for  the  small  (i.e.  kinetic  theory)  limit  whereas  fE,p  is  an  excess  function  of  importance  for 
larger  According  to  the  Stiel-Thodos  method. 


A, 


ref^p 


=  8.9775  X  10^ 


(Tc, 


V 


Zg,% 

Vc,p 


fE,p{Pred,p)  =  1-223  X  10  2[exp(0.535p,.grf  )  -  1]. 


(52) 

(53) 
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Properties  used  in  these  calculations  are  provided  in  Table  1  for  n-heptane,  iso-octane  and  all  light  species 
and  in  Table  2  for  the  quasi-steady  radicals.  For  these  latter,  estimates  of  pc  and  vc  are  made  using  the 
group  contribution  method  of  Joback  (Reid  et  al.*).  Then,  Tq  is  estimated  assuming  Zq  =  0.28;  the  same 
assumption  is  made  to  estimate  the  Pc,vc  and  values  in  parentheses  in  Table  1. 


III.  Numerical  method 


The  equations  are  solved  under  the  large  Peclet  number,  Pe,  assumption  which  is  consistent  with  using  the 
inviscid  momentum  equation.  Basically,  a  large  Pe  value  implies  that  convection  dominates  diffusion  which  is 
the  essence  of  the  problem  under  consideration.  If  one  defines  a  dimensionless  parameter  Si  =  —miJi/{pu)o, 
then  one  may  construct  =  —eijYi  =  pUD,i/ipu)o  where  UD,i  is  the  species  diffusion  velocity.  The  large 
Pe  value  means  that  |e'|  <C  1. 

Three  dimensionless  parameters  are  defined 


Tref  Aq  Ru 

along  with  characteristic  chemical  kinetic  rate-based  gradient  quantities 

niiTZi 


GYR,  = 
GTR  = 

and  a  characteristic  conduction  length 


{pu)o 
m 

Cp{pu)o  \ 


NcKh- 


1 


RuTref  ; 

•'  i£L 


It  = 


m\ 


Ru  (p^)o 

which  is  of  same  order  as  £i.  Then,  the  governing  equations  become 


ax  ax 


a  /  dT 


-A-l  -  -Ax  GTR  =  — 


d  AlrdT 


rdT  I  Gp^i  Gp^n 


m  \  dx 


dx  \  m  dx 


ax  ^ '  '  m  - 


i£L 


rtii  tuh 


(54) 

(55) 

(56) 

(57) 

(58) 

(59) 


The  idea  is  to  expand  the  terms  of  these  equations  in  function  of  Si  and  retain  only  the  lowest  order  terms. 
Doing  so  results  in 


dx 

=  Ax  {GYRi  +  GYDi) 

(60) 

df 

dx 

=  Ax  {GTR  +  GTD) 

(61) 

where  quantities  GYDi  and  GTD  are  listed  in  Appendix  B. 
The  equations  are  solved  using  a  stiff-equations  integrator. 


IV.  Results 

Preliminary  results  are  discussed  here,  obtained  for  A  =  1,  so  as  to  be  relatable  to  other  simulations 
described  in  the  literature  using  this  typically  made  assumption.  These  results  are  meant  to  probe  (1)  the 
impact  of  diffusion  with  respect  to  a  case  where  there  is  no  diffusion,  and  (2)  the  validity  of  the  large  Pe 
assumption  on  which  the  numerical  method  was  constructed.  All  computations  were  performed  for  a  fuel 
that  was  a  mixture  of  90%  iso-octane  and  10%  n-heptane,  with  {pu)o  =  0.8  g/cm^-s  and  stoichiometric 
equivalence  ratio,  i.e.  (j)  =  1. 

Figures  1  and  2  illustrate  the  results  without  diffusion  and  with  diffusion  and  Tables  4  and  5  both  pertain 
to  simulations  performed  including  diffusion.  Without  diffusion  (fig.  1),  the  calculation  proceeds  through 
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an  incubation  region  and  a  sharp  increase  starting  at  a  location  ,  Xfs,  defined  as  the  fiame  location  (for 
accounting  purposes,  Xfs  occurs  at  Tq  +  900  K).  After  that  location,  T  increases  within  a  short  distance 
to  almost  its  final  value  which  is  in  excess  of  2800  K.  Defining  Ufs  as  the  characteristic  velocity  at  Xfs, 
without  diffusion,  Ufs  =  120  cm/s  and  the  value  reached  for  the  asymptotic  T  value  is  250  cm/s.  With 
diffusion  (fig.  2)  it  was  verified  that  the  large  Pe  assumption  assumption  (i.e.  |e'|  ^  1)  holds  in  the 
incubation  region.  However,  as  the  flame  is  approached,  there  is  an  indication  that  this  assumption  is  no 
longer  valid.  This  finding  is  not  unexpected  since  the  fiame  is  inherently  a  small-scale  structure  involving 
very  large  gradients.  As  a  results  of  the  \e^\  <C  1  assumption  breakdown,  the  calculation  stalls  at  T  values 
or  X  positions  corresponding  to  past  the  fiame  initiation.  Values  of  Xfs  at  two  Tq  and  three  po  conditions 
are  listed  in  Table  4  and  the  corresponding  Ufs  is  listed  in  Table  5. 


V.  Conclusions 


A  model  has  been  developed  for  simulating  quasi-steady  quasi-one-dimensional  laminar  premixed  flames 
in  the  configuration  of  a  jet  injected  in  air.  The  model  has  been  derived  in  the  framework  consistent  with  it 
being  used  in  conjunction  with  a  reduced  chemical  kinetic  model  based  on  constituents  and  species.  As  such, 
the  formulation  includes  the  option  of  a  jet  variable-area  with  position,  a  complete  mass-diffusion  matrix, 
and  a  mixture  thermal  conductivity  computed  taking  into  account  the  global  constituent  and  all  species. 
Furthermore,  a  real-gas  equation  of  state  is  used. 

The  equations  are  solved  under  the  large  Peclet  number  assumption.  Results  show  that  this  assumption 
is  appropriate  in  the  flame  incubation  region  but,  not  unexpectedly,  it  breaks  down  very  near  the  fiame. 
Further  efforts  will  be  devoted  to  enlarging  the  numerical  technique  to  accommodate  both  the  incubation 
region  and  the  flame  region.  Simulations  performed  devoid  of  diffusion  show  that  the  numerical  method  can 
handle  the  strong  reaction  region  of  the  fiame  very  well. 
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Appendix  A 

Miscellaneous  relationships  relevant  to  the  EOS  are 


Q-n 


XpJCji  cipn  (T) 

5  ^mix  Xpbp 


p  n  p 

where  indices  do  not  follow  here  the  Einstein  notation,  and 

^pn  ~  (1  ^  )\/CX.ppCXnn', 

(Xpp{T)  =  0.457236(i?„Tc-^p)^  1  -I-  Cp(l  —  ^jT^ed,p)  IPc,p, 

Cp  =  0.37464  -h  1.542260^  -  0.269920^, 

where  Tred,p  =  T/Tc,p,  Tc,p  is  the  critical  temperature  and  Dp  is  the  acentric  factor.  Also, 

,RuTc,p 


bp  =  0.077796- 


PCd 


Tc,pn 


(1  kpn)  \/Tc^pPc,n  with  kpp  —  0, 

1  /  1/3  ,  l/3p 
VC,pn  =  g  +  Wc,n  j  > 

^C,pn  —  2  i^c,p  4"  ^C,n)  1 
Tc,pnZc,pn 

PC,pn  —  1 

'^C,pn 


(62) 

(63) 

(64) 

(65) 

(66) 

(67) 

(68) 

(69) 

(70) 
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with  Tred,pn  =  TjTc^pn,  Zc,p  being  the  critical  compression  factor  defined  as  Z  —  p/{pTRu/m),  vc,p  be¬ 
ing  the  critical  volume,  and  pc,p  being  the  critical  pressure,  kpn  is  an  empirical  mixing  parameter.  The 
relationship  between  the  parameters  kpn  and  is 


(1  kpn)  —  (1  kpn) 


{vc,pvc,n)^^‘^ 

'^C,pn 


(71) 


Values  of  kpn  are  listed  in  Table  3,  where  the  values  were  obtained  from  Reid  et  al.®  or  Knapp  et  al.^® 
Otherwise,  for  pairs  involving  O2,  N2,  CO,  CO2  or  H2O,  k^n  —  0. 


Appendix  B 

The  coefficients  of  eqs.  60  and  61  are 

GTD  =  X  GT R  X  (1  +  T  T 

GYD,  =  (73) 

jeLUNS 


where 


£t  = 


SR,i  = 


— 


Cp 

rmn 

{pu)o 

P 


+  miy] 

dx  ^  ^  5T 


rriH 


eR,r  +  {A)‘^  Y,  IcjxGYD, 

jeLUNS 


A 


TrefB^  i  X  GTR  +  Y  -  7?“) —  ^ 

jGLUNS 


Cp{pu)Q 


=  (A)2  X  GTR  X 


m 

rriH 

pT^ref 


rrij 


R  = 


Mij  _ 


dlnXl 


Cp{pu)o 

^£R,i 


y{g:^-—c-. 

\  niH 


p,H  I 


dx 
rriin  dA 
{pu)o  dx 


X  SR^i  +  Ax  GTR 


ST 


m 


(A,-21“)  — 


rrii 


In  eqs.  74-79,  the  following  replacements  were  made 


d  (lTdT\  ~5_{i^xGTR)dT 


dx  ym  dx 
dsRi 

by 


ST 


dx 


SsR^i  dT 


dx  ST  dx 

and  quantities  J  ()  /ST  are  numerically  calculated  using  finite  differences. 


(74) 

(75) 

(76) 

(77) 

(78) 

(79) 
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Species 

m  (g/mol) 

Tc  (K) 

PC  (bar) 

vc  (cm® /mol) 

n 

C7H16 

100.2 

540.2 

27.4 

432 

0.349 

CgHis 

114.23 

568.8 

24.9 

492 

0.396 

N2 

28.013 

126.26 

33.4 

89.8 

0.39 

H2O 

18.015 

647.3 

221 

57.1 

0.344 

CO2 

44.01 

304.1 

73.8 

93.9 

0.225 

O2 

32.0 

154.6 

50.43 

73.4 

0.025 

H 

1.008 

33.15 

(12) 

64.2 

(0.) 

H2 

2.016 

32.94 

12.84 

64.3 

-0.216 

CO 

28.01 

132.9 

35 

93.2 

0.066 

CH4 

16.04 

190.6 

45.94 

99.2 

0.0108 

H2O2 

34.015 

728 

220 

(77) 

(0.1) 

C2H2 

26.04 

308.3 

61.4 

112.7 

0.19 

C2H4 

28.054 

282.4 

50.4 

130.4 

0.089 

CH2O 

30.026 

408 

65.9 

(144) 

0.253 

Table  1.  Fuel  and  unsteady  light  species  properties. 


Radicals 

m  (g/mol) 

Tc  (K) 

PC  (bar) 

Vc  (cm® /mol) 

n 

0 

16.0 

116 

76 

35 

0. 

CH 

13.02 

183 

73 

58 

0. 

CH2 

14.027 

210 

66 

73 

0. 

CH3 

15.034 

221 

62 

82 

0. 

OH 

17.01 

167 

85 

45 

0. 

HCO 

29.02 

299 

70 

99 

0.25 

OOH 

33.01 

226 

83 

63 

0. 

HC2 

25.03 

291 

67 

100 

0.2 

C2H3 

27.044 

292 

57 

119 

0.1 

Table  2.  Quasi-steady  radicals  properties. 
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p 

n 

k' 

alkane 

alkane 

0.0 

alkane 

N2,  O2 

0.15 

alkane 

CO2 

0.11 

alkane 

H2O 

0. 093-0. OOGhc 

alkane 

H2 

0.099nc-0.81 

H2 

N2,  O2 

0.12 

H2 

CO2 

-0.162 

CO2 

H2O 

0.095 

CO2 

N2,  O2 

-0.017 

H2O 

N2,  O2 

0.17 

Table  3.  Values  of  k  for  species  pairs,  ric  is  the  number  of  C  atoms  in  the  species. 


Po,  bar 

20 

30 

40 

T,  K 

700 

2.39 

1.36 

0.9 

800 

0.9 

0.4 

0.2 

Table  4.  Flame  position,  Xfs  (cm)  for  stoichiometric  conditions  and  simulations  including  diffusion. 


Po,  bar 

20 

30 

40 

T,  K 

700 

118 

88 

64.4 

800 

132 

91 

73.4 

Table  5.  Characteristic  speed  (cm/s)  at  Xfs  for  stoichiometric  conditions  and  simulations  including  diffusion. 
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(jUi/ioiu)  “n  (jUI/ioui)  '°n  (jUi/IOUi) 


Figure 


. . 

0  0.1  0.2  0.3  0.4  0.5  0.6 

b  *  (cm) 


.  Predictions  of  the  reduced  model  without  diffusion  at  pQ  =  30  bar,  Tq  =  800  K  and  (p  = 
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(mol/nn“)  N^2  (mol/nn“)  (mol/m’) 


. . 

0.1  0.2  0.3  0.4 

X  (cm) 


oooooooooooo^ 


Pq  =  30  bar 
T.  =  800  K 
(b=1 


0.1 


o 

o 


o 

o 


o 

o 

© 

o 

o 

o 

o 


0.2 

X  (cm) 


0.3 


0.4 


10 

9 

8 

7 

6 

5 

4 


Pq  =  30  bar 
T.  =  800  K 
(^=1 


O 

O 


0.04 


Pq  =  30  bar 
T.  =  800  K 
(^=1 


O 

o 

o 

o 


0^  oo  oc^^oo  ooj^QOcoGay^a! 

d  X  (cm) 

2.5r 


2h 


p„  =  30  bar 
Tq  =  800  K 

(1)  =  1 


DOOOOOOOOOOOOO  'telSSSSCOOO  o  o 
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Modeling  of  Steady  Laminar  Flames  for 
One-dimensional  Premixed  Jets  of  Heptane/Air  and 

Octane/ Air  Mixtures 

Kenneth  G.  Harstad^’*and  Josette  Bellan^’**d 
Jet  Propulsion  Laboratory* ,  California  Institute  of  Technology, 

Pasadena  CA  91109-8099 

California  Institute  of  Technology** ,  Mechanical  Engineering,  Pasadena,  CA  91125 


A  model  is  proposed  for  quasi-one-dimensional  steady  flame  development  in  the  con¬ 
figuration  of  an  inviscid,  premixed  fuel  jet  injected  into  air.  The  governing  equations 
are  written  within  the  framework  of  a  reduced  kinetic  model  based  on  constituents  and 
species.  The  reduced  kinetic  model,  previously  exercised  in  a  constant-volume  perfectly- 
stirred  reactor  mode,  has  been  successful  at  predicting  ignition  and  combustion  product 
and  temperature  evolution  for  n-heptane,  iso-octane,  PRF  fuel  combinations,  and  mixtures 
of  iso-octane  with  either  n-pentane  or  iso-hexane.  The  differential  governing  equations  have 
the  option  of  an  axially  variable  area  and  they  are  coupled  with  a  real  gas  equation  of  state. 
The  flame  development  model  accounts  for  a  full  diffusion  matrix,  and  thermal  conductiv¬ 
ity  computed  for  the  species  mixture.  Results  from  four  simulations  at  various  conditions 
are  presented. 


I.  Introduction 

To  predict  the  initiation  and  evolution  of  turbulent  reacting  flows  in  practical  devices  one  must  couple 
chemical  kinetics  and  turbulence  models.  Because  codes  resulting  from  these  models  are  very  computationally 
intensive,  the  chemical  kinetic  mechanisms  cannot  be  utilized  in  their  entire  complexity  and  instead  they  are 
reduced  to  mechanisms  that  are  sought  to  be  compact  as  well  as  reliable.  Recently  a  reduced  kinetic  model 
relying  on  constituents  and  species^^  has  been  developed,  which  was  successful  in  reproducing  ignition  time, 
and  combustion  product  and  temperature  evolution  for  n-heptane,  iso-octane,  PRF  fuel  combinations,  and 
mixtures  of  iso-octane  with  either  n-pentane  or  iso-hexane. 

In  this  conceptual  model,^’^  the  ensemble  of  species  in  a  full  (or  skeletal)  mechanism  is  partitioned  into 
heavies  (carbon  number,  n  ^  3)  and  lights  (the  remaining  of  the  set).  The  heavies  can  be  either  radicals  or 
stable  species.  The  lights  are  oxygen,  nitrogen,  the  final  combustion  products  and  light  radicals/molecules 
(e.g.  CH3,  CH4,  H2O2).  The  heavies  are  not  treated  as  a  species  set,  but  as  a  set  of  base  constituent  radicals 
which  correspond  to  radicals  as  used  in  group  additivity  theory.^  Sometimes  light  species  are  radicals  which 
may  have  the  same  chemical  formula  as  constituents,  however,  the  difference  between  constituents  and  these 
light  species  is  that  the  later  are  unbound  to  other  chemical  entities  whereas  the  former  are  bound  to  other 
chemical  entities,  i.e.  other  constituents.  For  each  constituent  k,  its  molar  density,  Nk,  is  the  sum,  over  all 
heavy  species,  of  the  count  of  the  constituent  in  each  heavy  species  multiplied  by  the  molar  density  of  that 
species.  The  constituents  are  not  just  based  on  atom  counts  and  their  element  composition  is  not  linearly 
independent,  however,  the  constituents  are  linearly  independent.  Thus,  the  constituents  are  independent 
structural  elements  which  have  individual  valence  bond  topologies. 

For  each  heavy  species,  and  thus  for  the  entire  set  of  heavy  species,  there  is  a  unique  and  complete  set  of 
constituents.  The  subset  of  constituents  that  have  a  very  small  contribution  to  the  total  constituent  count 
having  molar  density  N^,  is  replaced  by  the  complementary  subset  of  constituents  already  accounted  for  in 

^Senior  Engineer. 

'^’Senior  Research  Scientist,  AIAA  Fellow  (corresponding  author,  josette.bellan@jpl.nasa.gov). 
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Nc-  The  rule  for  replacement  is  that  each  constituent  having  a  small  contribution  to  Nc  is  replaced  by  one 
having,  first,  close  valence  structure  and  second,  close  composition  to  the  constituent  neglected.  After  this 
replacement,  one  obtains  an  optimal  Nk  set.  It  is  this  optimal  set  of  constituents  that  is  used  in  the  kinetic 
reduction. 

To  test  the  capability  of  the  reduced  kinetic  model  in  the  a  simple  flow  model,  the  first  step  is  to  evaluate 
it  in  a  steady,  quasi-one-dimensional  configuration  before  undertaking  the  more  arduous  task  of  the  counter¬ 
flow  configuration  for  which  the  overwhelming  experimental  data  has  been  obtained.  We  first  present  the 
governing  equations,  as  they  are  more  general  than  the  typical  ones  for  steady  quasi-one-dimensional  flow; 
in  particular  we  use  here  accurate  transport  properties  and  a  real-gas  equation  of  state.  Further,  we  discuss 
some  results  obtained  with  the  model  and  then  present  a  concise  summary  and  conclusions. 


II.  Governing  equations 

The  configuration  of  interest  is  that  of  a  steady,  round  jet  in  a  system  of  coordinates  {x,  r,  0)  where  it 
is  assumed  that  there  is  uniformity  in  the  0  direction  and  thus  the  coordinates  of  interest  are  {x,r).  The 
entrance  of  the  jet  {x  =  0)  is  labeled  by  the  subscript  0  and  the  far  field  boundary  is  at  a:  =  L.  Within  the 
quasi-one-dimensional  framework,  all  quantities  are  averaged  in  the  r  direction,  and  thus  the  jet  boundary 
is  located  at  r  =  R{x).  Therefore,  the  jet  area  is  A{x)  =  'xR^{x).  A  quantity  ^{x)  =  Ro/R{x)  is  defined  to 
track  the  jet  spacewise  variation. 


A.  Continuity  and  momentum 


We  define  the  mass  density  p,  and  the  axial  and  radial  velocities,  u  and  v.  Radial  averages  are  designated 
by  0,  as  for  example 


p{x) 


R 

^  j  P{x,r)  r  dr. 
0 


(1) 


In  the  entire  derivation  below,  all  quantities  will  be  considered  as  radial  averages,  but  the  ()  symbol  will  be 
omitted  for  notation  simplicity. 


1.  Continuity  equation 

•  • 

Defining  the  mass  flow  rate  M,  the  steady  flow  assumption  implies  puA  =  M  =  constant,  which  means  that 


pu  =  M/A  =  M/wR^ix). 


Replacing  eq.  2  in  the  steady  continuity  equation 


1  djrpv)  d{pu)  ^ 

r  dr  dx 

leads  to 

M  dR 


(2) 

(3) 

(4) 


2.  Momentum  equation 

Under  the  assumption  of  an  inviscid  fluid,  the  steady  momentum  equation  components  are 


f  d  d\ 

/a  a  \ 


dp 

dx 

dp 


0 

0. 


(5) 

(6) 


To  make  progress  in  obtaining  an  equation  for  the  jet  area,  we  expand  p  and  p  ^  in  even  powers  of  r  (due  to 
symmetry  conditions),  as  a  means  of  taking  advantage  of  the  configuration  symmetry.  Thus,  we  postulate 
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that 


^  '^max  ^2n 

P  =  PG  +  l^Poul  ^n{x)^, 

n=0 


n=0 


2n 


where  cr„  and  /„  are  arbitrary  functions.  Replacing  these  expressions  in  eqs.  5  and  6  yields 

1  d(7n  _  ^2d{ffu)  ,  d  ^ 

2  dx  dx  ^  dx  dx  1 


dx 


di 


nan  =  RQ^—[fn-i—]  for  n  >  0. 


dx 


Combining  eqs.  9  and  10  for  n  =  1,  leads  to 

,d{efi) 


2r 


+  Rl 


dx  '  ^  ^  dx 


fo 


d^ 


dx\  dx 


=  0. 


(7) 

(8) 


(9) 

(10) 


(11) 


An  analysis  of  eq.  11  provides  insights  from  which  further  simplifications  are  obtained.  For  example,  near  a 
stagnation  point,  which  is  obtained  for  $  0,  the  second  term  on  the  left  hand  side  of  eq.  11  dominates  the 

first  term  which  can  then  be  neglected.  Also,  near  the  a:  =  0  boundary,  the  jet  is  nearly  radially  uniform, 
which  means  that  /i  =  0.  These  simple  considerations  indicate  that  /i  can  be  neglected  over  the  entire 
length  of  the  jet.  To  analyze  the  last  term  of  eq.  11,  we  define 

=  I  —  k  [  (^{x')— — ^-x'dx'  (12) 

Jo  Po 

where  k  is  determined  by  the  downstream  boundary  condition  and  is  therefore  constant.  Then,  the  last  term 
of  eq.  11  contains 


f  Po!^\ 

dx  \  p  dx  J 


k  2  dC  k 
'2^^dx  2^ 


l  +  k( lx\-^  -  r  ax')^x'dx' 

V2  Po  Jo  Po 


(13) 


Outside  of  the  flame,  {p/po)  ~  1)  which  together  with  the  fact  that  the  right  hand  side  of  eq.  11  must  be 
constant,  leads  to  assuming  that  (  is  constant.  This  implies  that  the  right  hand  side  of  eq.  13,  i.e.  {—k(^/2) 
is  also  nearly  constant.  To  remove  the  ambiguity  of  dealing  with  the  product  of  two  constants,  we  take  C  =  1 
which  is  a  good  approximation  for  satisfying  eq.  11.  In  the  flame  region,  p/pQ  <  1;  for  fixed  p,  its  typical 
(or  mean)  value  is  that  of  To/T  «  1/4  found  in  flames.  However,  the  flame  thickness  is  small,  and  a  value  of 
C  ^  1  but  varying  weakly  may  be  used  in  this  region  since  the  right  hand  side  of  eq.  11  will  still  be  nearly 
constant.  For  the  results  presented  in  section  IV,  the  deviation  of  C  from  unity  was  0(10^^). 

For  A  constant,  none  of  these  considerations  apply  and  the  momentum  equation  becomes  a  vehicle  to 
compute  Vp,  but  under  the  assumption  of  very  small  Mach  number  the  pressure  is  assumed  constant  for  the 
calculation  of  the  equation  of  state  (EOS). 


B.  Species  and  energy  equations 

The  complete  formulation  for  the  species  and  the  energy  equations  has  been  presented  elsewhere."*  However, 
that  formulation  was  for  species  and  it  is  not  applicable  here  where  the  model  is  for  both  constituents  and 
species.  In  this  reduced  kinetic  model,  the  evolution  of  the  global  constituents  molar  density  is  found  from 
input  tables  and  it  is  only  the  light-species  equations  that  are  progress  variables,  together  with  the  energy; 
however,  the  constituents  make  contributions  to  the  light,  and  these  contributions  must  be  now  formulated 
in  the  species  equations  in  terms  of  the  constituents.  To  this  end,  we  assign  indices  i,j  and  q  to  an  ensemble 
L  denoting  the  lights,  k  to  the  ensemble  C  of  the  constituents  and  l,m  to  the  ensemble  H  of  the  heavies. 
Ensemble  L  is  further  partitioned  into  that  of  the  quasi-steady  species,  LQS,  and  that  of  the  unsteady 
species,  RUNS.  All  quantities  are  assumed  averaged  over  the  r  direction,  so  there  is  only  a  variation  with  x. 
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1.  Species  equation 

We  partition  the  species  molar  flux  in  the  light-species  governing  equation  in  two  categories  as  Ji^ught  which 
represents  the  contributions  on  light  species  i  of  all  lights,  and  Ji^teavy  which  represents  contributions  on 
light  species  i  from  all  heavies 


Ji  —  Ji, light  Ji, heavy 

(14) 

/  __  _  d\nT  dXj\ 

+Y^-^  dx  j 

(15) 

V  "T  dXl 

heavy  —  ^  ^  J^il  —  JH,i 

(16) 

where  n  is  the  molar  density,  X  is  the  mole  fraction,  T  is  the  temperature  and  quantities  Vij  are  pairwise 
diffusion  coefficients  for  the  lights  that  are  computable  from  the  binary  diffusion  coefficients  through  mixing 
rules;  ways  to  compute  binary  diffusion  coefficients  have  been  given  by  Harstad  and  Bellan®  and  mixing  rules 
were  derived  by  Harstad  and  Bellan.®  The  term  proportional  to  VT  includes  the  entire  Soret  effect  (in  Ji 
the  term  proportional  to  Vp  is  considered  negligible  because  it  is  usually  small  and  we  additionally  make  the 
small  Mach  number  approximation)  and  the  Soret  effect  is  included  in  Jijight  because,  as  for  'Dij,  Dx^i  can  be 
computed  using  mixing  rules.  These  mixing  rules  generally  include  contributions  from  the  heavies,  however, 
since  the  individual  heavies  are  not  individually  tracked,  we  use  a  single  mean  heavy  species  dynamically 
computed  from  the  injected  alkane  fuel  as  it  pyrolyses.  The  goal  is  now  to  express  JH,i  as  a  function  of  a 
constituent  equivalent  of  a  mole  fraction. 

If  one  considers  a  constituent  Nk,  then 

Nk  =  J2  N,  =  J2Nk  (17) 

l  k 

where  C^i  is  the  count  of  constituent  k  in  species  1.  To  find  a  way  of  writing  the  heavies  as  a  function  of  the 
constituents,  one  must  invert  the  matrix  of  elements  Cki-  However,  because  there  are  fewer  equations  (i.e. 
constituents)  than  unknowns  (i.e.  heavies),  the  inversion  can  only  be  approximately  performed.  To  perform 
this  inversion,  we  select  the  Householder  transformation  which  when  utilized,  yields 

^  ^  Cmk^k  where  6>fYii  =  ^  ^  (1^) 

k  k 


with  Smi  being  the  Kronecker  symbol.  One  may  interpret  Nm  as  the  minimum  norm  of  heavy  molar  densities 
for  given  Nk-  This  norm  favors  heaviest  species  and  it  is  thus  most  accurate  at  early  times,  when  heavy  species 
diffusion  is  important.  Since  Xi  —  Ni/n, 


JH,i  =  —  ^  ^ T’iiCik 


dNk 

dx 


Y^VuNi 


dln(n) 

dx 


(19) 


Since  ultimately  only  W  =  Nk  is  available  in  the  model,  an  additional  transformation  is  needed  in  the 
expression  of  eq.  19.  Under  the  verified  assumption  that  individual  dominant  constituents  of  molar  density 
Nk  are  mostly  quasi-steady,  it  was  postulated^  that  there  exists  a  mole  fraction  XCk  such  that  Nk  —  NcXCk- 
Therefore, 


XCk  = 

Lk  I 


(20) 


Finally,  it  is  necessary  to  compute,  from  the  "Du  coefficients,  a  global  diffusion  coefficient  for  the  term 
accounting  in  the  formulation  for  the  global  constituent  molar  density.  For  this  purpose,  we  define 


T^H,i 


1 

n 


Y'l’M  =Y'^aXi  =  XhV,h 

I  I 


(21) 
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where  Xh  is  the  global  mole  fraction  of  the  heavies  and  ViH  is  a  mean  diffusion  coefficient  from  the  light 
species  i  to  the  ensemble  of  the  heavies.  If  one  defines  a  mean  total  constituent  count  as 

n 


Cave=  —  ^y^CklXu 
k,l 


(22) 


then  one  may  approximate  JH,i  by 


JH,i  =  —nDn, 


This  leads  to 


Ji  =  -n 


dlnCg 

dx 

dlnCg. 
dx 


-  XiDta 


dluT 

dx 


(23) 


(24) 


To  write  Ji  as  a  function  of  mass  fraction  gradients  and  take  advantage  of  the  partition  of  the  lights  set  into 
unsteady  and  quasi-steady  species,  we  now  assign  indices  i  and  j  to  the  unsteady  light  species  and  q  to  the 
quasi-steady  light  species.  Thus  the  total  unsteady  mole  fraction  is 


Xg  =  l-XH-Y,X, 


and  the  mixture  mean  mass  becomes 


Xg 


m  — 


Yi 


(25) 

(26) 


Introducing  these  definitions  in  eq.  24  yields 
Jr  = 


~n 


=  —n 


dx 


dXg  ^  dlnC, 

/)  TJtq  +  XhTJiH 


dx 


dx 


X^Dt, 


-J^a  dXg 

*  dx 


\ ^  \  ^  dXi 

^(A,-2)“) - ^ 


rrij  dx 


d\nT 
dx 

9  I  Y  T"  dlnCg 
dx  dx 


dXg 


■  XiDT,i 


d\nT 

dx 


(27) 


having  used  dXj  =  {m/mj)dYj  -I-  Xjd\n{m)  and  having  defined  a  mean  unsteady-lights  mass  diffusion 
coefficient 

(28) 

^  j 

Writing  eq.  27  consistent  with  the  kinetic  reduction  modek’^  in  which  variables  Xq^Xn  and  Cggg  are  ob¬ 
tained  from  the  reduced  model  as  functions  of  the  temperature,  so  that,  for  example  dX^f/da;  =  {5Xfi/5T){dT/dx), 
leads  to 


where 


Ji  =  -n 

Dt 


SXg 


dx 

din  a. 
5T 


■n 


gSXn 

ST 


(29) 

(30) 


acts  as  an  effective  Soret  coefficient  for  the  reduced  model. 

Denoting  the  total  reaction  rate  of  species  i  hy  TZi,  the  governing  equation  for  species  of  mass  fraction 
Yi  is 


dY,  A  d  A  ' 

=  —m,n,  -  —  —m,J, 
ax  Vm  / 


In  accordance  with  the  reduced  kinetic  model 

TZi  —  TZ 


I  heavies 


(31) 


(32) 


lights  dZi  I 

where  TZi\ngi^^^  is  computed  using  CHEMKIN  II  in  conjunction  with  the  LLNL  rate  data,  and  | /lea-uies  ~ 
NgKnety  where  Knety  is  computed  by  interpolation  from  tables  where  it  is  listed  as  a  function  of  T  for  fixed 
(To,poj<(')  with  4>  being  the  equivalence  ratio  and  subscript  0  denoting  initial  conditions. 
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2.  Energy  equation 


We  write  the  energy  equation  under  the  assumption  that  the  Dufour  effect  has  a  negligible  contribution  to 
the  heat  flux.  Additionally,  since  in  our  model  the  heavy  species  are  combined  into  a  single  constituent, 
we  approximate  the  molar  enthalpy  of  all  heavy  species  averaged  over  the  molar  fluxes  by  the  enthalpy 
hn  =  {1/Xh)  Xli  which  is  the  mean  heavy-species  enthalpy  based  on  an  average  of  heavy  species  molar 
enthalpies  hi.  The  result  is 


E 


mj  ■ 
mH 


‘if 


J^ 


(33) 


where  A  is  the  mixture  thermal  conductivity.  The  energy  equation  is  thus 


CpdT  _  _d 
m  dx 


Cp^H  \ 
rriH  ) 


{rriiJi) 


+  'y  ^  hilZi  —  RuTrefNcK 

i 


(34) 


where  Cp  is  the  heat  capacity  at  constant  pressure  and  Kh  is  computed  by  interpolation  from  tables^  where 
it  is  listed  as  a  function  of  T  for  fixed  (To,po,(j)). 


C.  Equation  of  state 

The  pressure  is  calculated  from  the  Peng-Robinson  (PR)  EOS 

_  RuT  a^jiix 

{vpR  -  b  mix^  (^Pi?  +  ‘ibmixVpR 

where  vpR  is  the  molar  volume  and  Omix  and  hmix  are  functions  of  T  and  W  (see  Appendix  A). 


(35) 


D.  Transport  properties 

To  derive  transport  properties  consistent  with  the  kinetic  model  approach,  we  define  an  ensemble  of  all 
light  species  and  one  mole-averaged  heavy  species:  LU  H.  Indices  p,  n  and  r  refer  to  any  species  from 
this  ensemble.  The  transport  properties  under  consideration  are  the  diffusion  coefficients  Rpn,  the  thermal 
diffusion  factors  Dpp,  and  the  thermal  conductivity;  the  viscosity  does  not  enter  the  calculations,  since  an 
inviscid  situation  is  assumed,  however,  because  transport  property  calculations  of  thermal  conductivity  and 
viscosity  are  very  much  related,  we  borrow  from  methods  to  compute  viscosity  in  order  to  compute  the 
thermal  conductivity. 

For  the  diffusion  coefficients,  the  first  task  is  to  compute  the  binary  diffusion  coefficients  which  are  the 
building  blocks  of  the  pairwise  diffusion  coefficients.  To  this  end,  we  adopt  the  method  of  Harstad  and 
Bellan®  which  gives  (in  cgs  units) 


77  ]~) 
itj-ypn 


=  2.81  X  10'® 


fpiT) 

2/3 


(36) 


where  the  subscript  C  denotes  the  critical  state,  /d(T)  =  {TredY  with  Ins  =  X)a=o  where  the 

a®  vector  has  elements  {-0.84211,  -0.32643,  -0.10053,  0.07747,  0.0127,  -0.00995},  and  rp  is  a  constant  0(1) 
which  provides  an  empirical  adjustment  for  the  specifics  of  the  collisional  interactions  of  a  selected  pair  of 
species.  Values  of  rp  are  listed  in  Harstad  and  Bellan®  for  species  pairs  relevant  to  combustion.  The  heavy 
species  molar  mass  is  part  of  the  reduced  model  tables  and  decreases  as  the  reaction  proceeds. 

From  the  Dpn  ’s,  the  pairwise  diffusion  coefficients®  can  now  be  computed  using,  as  an  approximation  (a 


6  of  14 


American  Institute  of  Aeronautics  and  Astronautics 


truncated  series  proposed  by  Ern  and  Giovangigli^), 


V  — 


Dr>n  — 


VI  = 


Xp  ^  ^  DpnCXjj 

,nr ; 
n 

(i  +  ^p) 


'DTD, 


N 


Xr, 


D;5p^  +  (1  -  +  '^nD:)  +  J2iYrarD:), 


,  -‘-^pn 


=  'T(iyy.(+Yy. 


N 


n=tp 


D* 

Drti 


where  the  mass  diffusion  factors  are  computed  as 


_  dXp  dlii-ip 

O-D.pn  —  Qj^  + 


dX^ 


(37) 

(38) 

(39) 

(40) 

(41) 


from  the  EOS,  with  7^  =  (ppl(p°,  where  (f  is  the  fugacity  coefficient  and  the  superscript  o  denotes  the  pure 
(Xp  =  1)  limit. 

According  to  Harstad  and  Bellan® 


Dt,p 


to 


T 

P 


N  /  N 


E  E^^c.  n{Tred,r  n) 


n—1  \  r=l 
\r^n 

(^p  ^p 

R^uTi 


{mrCv'^  -  TTlnC^J) 
{rrir  +  mn)Dr  - 


^nVnpf 


(42) 

(43) 


where  Xp  is  the  species  p  thermal  conductivity  computed  according  to  eq.  51  and  ui^  is  a  weighting  factor 
given  by  eq.  46.  Factors  Crra(^red,rn)  are  closely  related  to  the  temperature  derivative  of  Dp„  taken  at 
constant  pressure,®  and  they  are  computed  from 


Cr  n(pred,r  n)  —  ^  g 


Sr  n  +  Tred,  r  n  ffi(7r.ef;^7.  T) 


dSr 


dTr 


red,r  n  J 


where  Tred,r  n  is  defined  in  Appendix  A. 

To  compute  A,  we  use  the  Wassiljewa-Mason-Saxona  method  described  in  Reid  et  al.,^ 


(44) 


X„c^QA„  (45) 

n^LVJH 

where  the  weighting  factors  are  computed  as  recommended  in  Reid  et  ah,® 


-1 


</>pn  = 


(46) 


(47) 


where  rjp  represents  the  viscosity  of  species  p.  To  compute  the  individual  species  viscosities,  we  adopt  a 
method  explained  in  Reid  et  al.®  whereby 


Vp  Vref,p  f R—TiYred.pj  1 


(48) 


7  of  14 


American  Institute  of  Aeronautics  and  Astronautics 


where  the  subscript  R  —  T  stands  for  “Roy-Thodos”  and 


fR-T{Tred,p)  =  2.25(exp(0.0464Tred,p)  -  exp(-0.2412rred,p)), 

/  z  \ 

=  1.08x10^4(  )i/2  ^  ^ 

V  ^c.p  / 


(49) 

(50) 


The  other  ingredient  entering  the  A  calculation  is  the  expression  for  Ap  which  is  here  computed  (in  cgs  units) 
using  the  Stiel-Thodos  method  (Reid  et  al.®)  as 


~  ^ref,p  \f R~T(Tred,p)  +  ^5  f E,p{Pred,p)^  ) 


/a  =  3.75 


R^ 


(51) 


(52) 


0.7862  -  0.71090  +  1.316802 

where  f\  is  a  fnnction  of  the  acentric  factor  0  given  by  the  Chung  at  al.  formula,®  fR^T{Tred,p)  of  eq.  51 
accounts  for  the  small  (i.e.  kinetic  theory)  limit  whereas  Je^p  is  an  excess  function  of  importance  for 
larger  According  to  the  Stiel-Thodos  method, 


A, 


■ref,p 


=  8.9775  X  10^ 


fTcd 


V  mp 


1/2 


Zc,p 

Vc,p 


2/3 


fE,p{Pred,p)  =  1-223  X  10  2[exp(0.535ppgrf_p)  -  1]. 


(53) 

(54) 


Properties  used  in  these  calculations  are  provided  in  Table  1  for  n-heptane,  iso-octane  and  all  light  species 
and  in  Table  2  for  the  quasi-steady  radicals.  For  these  latter,  estimates  of  pc  and  vc  are  made  using  the 
group  contribution  method  of  Joback  (Reid  et  al.®).  Then,  Tq  is  estimated  assuming  Zq  =  0.28;  the  same 
assnmption  is  made  to  estimate  the  pc^vc  and  values  in  parentheses  in  Table  1. 

III.  Numerical  method 

The  equations  are  solved  using  a  numerical  method  relying  on  three  to  fonr  steps,  as  follows.  In  a  first 
step,  Based  on  prohles  of  T  and  Y^’s  obtained  from  published  experimental  data  and  models  available  in  the 
literature,  an  initial  guess  is  made  for  these  profiles  as  a  function  of  x,  for  a  specified  grid.  In  the  second  step, 
a  preliminary  relaxation  integration  relative  to  the  grid  is  made  by  combining  the  assumed  profiles  of  the  first 
step  with  that  obtained  using  differences  from  the  conservation  eqnations.  Results  are  examined  to  assess 
if  this  preliminary  relaxation  procednre  gives  reasonable-looking  profiles;  if  not,  the  first  step  is  revised.  In 
a  final  step,  further  profile  relaxations  are  made  using  a  pseudo-unsteady  time  increment  procedure.  This 
procedure  is  meant  to  relax  the  numerical  residues  of  the  conservation  equations.  For  variables  Yi  and  T, 
the  increments  are 


AF, 

AT 


where 


uAt 

uAt 


Cpe  =  Cp 


'rriiTZi 

--I 

(  m^JA 

_ 

pu 

dx 

[  pu  ) 

dx 

■  d 

(  A  dT\ 

1 

m 

dTl 

dx 

\pu  dx  / 

pu 

Cpe 

dx  j 

m 


pu 


r  -—r 

^pi  ^pH 

rriH 


A. 


(55) 

(56) 

(57) 


The  axial  coordinate  is  discretized  in  intervals  Ax  =  0(10^^)  cm  and  the  iteration  pseudo  steps  are  At  = 
0(0.1)  ps.  Since  u  —  O(IO^)  cm/s,  each  time  iteration  pseudo  step  corresponds  to  an  axial  distance  change 
of  O(10~®)  cm,  or  approximately  1%  of  Ax.  For  high-p  flows,  a  fourth  step  is  necessary  because  of  the  very 
steep  flame:  we  re-grid  successively  in  the  flame  region  in  order  to  obtain  good  resolution. 

IV.  Results 

Several  sets  of  results  are  discussed  here,  some  of  which  are  comparable  with  published  experimental 
data.  Because  experimental  data  is  not  usually  accompanied  by  illustrations  of  initial  condition  profiles,  and 
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because  the  computations  tend  to  be  sensitive  on  how  close  the  initial  profiles  are  to  the  quasi-steady  flame 
ones,  simulations  are  still  challenging  from  the  numerical  viewpoint  in  terms  of  achieving  profile  convergence 
under  the  one-dimensional  assumption. 

A.  Case  1 

The  first  case  is  one  for  which  experimental  data  exists.^®  Experiments  were  performed  for  a  premixed 
n-heptane/air  jet  under  the  following  conditions:  (j)  =  1,  p  —  1  bar,  Tq  =  350  K,  L  =  0.6  cm,  {pu)o  —  0.18 
g/cm^-s.^°  Illustrated  in  Fig.  1  are  computations  using  our  reduced  kinetics  model  (Fig.  la)  as  well  as 
experimental  results  (Fig.  lb)  .  Clearly,  the  flame  location  is  reproduced  with  good  accuracy,  the  velocity 
up  to  the  location  of  (du/dx)  =  0  is  very  well  predicted  but  its  computed  value  after  the  flame  exceeds  that 
of  the  experimental  data;  the  temperature  is  also  very  well  predicted  up  to  past  the  flame,  however,  the 
value  at  which  levelling-off  occurs  is  higher  than  in  the  experiment.  Given  the  compact  reduced  kinetics, 
some  inaccuracies  are  unavoidable.  Further,  in  Fig.  Ic  are  the  reactant’s  profiles.  The  mass  fractions  of  the 
reactants  are  constant  until  the  flame  region  where  a  reduction  in  1^02  is  observed  concomitant  with  a  small 
increase  in  the  Yh;  the  Fo2  reduction  is  attributed  to  the  initiation  of  the  reaction  whereas  the  small  increase 
in  Yh  is  attributed  to  diffusion  promoted  by  the  large  temperature.  The  same  diffusion  process  explains 
the  further  increase  in  Yo2  until  the  kinetic  rates  take  over  and  both  Yh  and  y’o2  decrease  precipitously. 
Noteworthy,  the  heavy  species  disappear  before  all  oxygen  is  exhausted,  and  the  remaining  oxygen  is  used 
for  the  light  species  reactions,  in  particular  to  convert  CO  to  CO2. 

B.  Case  2 

The  simulations  of  Case  2  are  devoted  to  an  iso-octane/air  mixture  under  the  conditions  of  </>  =  1,  p  =  10  bar. 
To  =  300  K,  T  =  0.4  cm,  (pu)o  =  3.0  g/cm^-s.  Results  are  displayed  in  Fig.  2.  A  quantitative  comparison 
with  the  n-heptane/air  flame  must  take  into  account  the  combined  effect  of  the  lower  initial  temperature  but 
larger  initial  pressure,  shorter  domain  and  larger  flow  rate.  Both  velocity  and  temperature  are  larger  than 
in  Case  1,  and  also  the  first  decrease  in  y/2  and  increase  in  Y^  are  muted. 

C.  Case  3 

For  Case  3,  we  are  returning  to  n-heptane/air  mixtures  but  at  a  different  initial  conditions  than  Casel.  For 
Case  3,  (p  =  1,  p  —  1  bar.  To  =  500  K,  T  =  0.6  cm,  (pu)o  =  0.10  g/cm^-s  and  the  results  are  exhibited  in 
Fig.  3.  By  comparison  with  Case  1,  To  is  here  larger  but  (pu)o  is  smaller.  Expectably,  the  temperature  is 
larger,  and  it  is  conjectured  that  this  larger  temperature  induces  the  larger  velocity.  The  reactants’  profiles 
are  qualitatively  similar  to  those  of  Case  1,  and  the  persistence  of  oxygen  farther  downstream  than  in  Case 
1  is  noted. 

D.  Case  4 

To  show  the  capabilities  of  the  model  for  rich  flames,  in  Case  4  we  compute  a  flame  similar  but  not  identical 
to  that  which  has  been  experimentally  addressed  by  El  Bakali  et  al.^^  The  reason  that  the  same  flame  could 
not  be  computed  is  that  for  the  very  small  mass  flow  rate  in  the  experiment,  the  code  did  not  converge; 
moreover,  the  domain  size  used  in  the  experiment  is  unclear  except  for  the  fact  that  L  >  0.4  cm.  Thus, 
the  initial  conditions  of  the  computation  are:  </>  =  1.9,  p  —  1  bar,  Tq  =  300  K,  L  —  0.5  cm,  (pu)o  =  0.0147 
g/cm^-s  (the  mass  flow  rate  is  here  increased  by  a  factor  of  approximately  2.5  from  that  in  the  experiment). 
Results  are  depicted  in  Fig.  4.  Both  the  temperature  and  the  velocity  have  smaller  values  than  those  of 
Case  1  where  Tq  was  slightly  larger,  the  mass  flow  rate  was  larger  by  an  order  of  magnitude  than  in  Case 
4,  and  the  domain  size  was  slightly  larger,  but  the  general  shape  of  the  profiles  is  qualitatively  similar  in  all 
cases.  Similar  to  the  experiment, the  results  of  Fig.  4  show  substantial  T/o  >  Yco2  (the  experimental  data 
shows  mole  fractions,  whereas  we  display  mass  fractions),  and  y/20  >  Yco2- 

V.  Summary  and  conclusions 

A  model  has  been  developed  for  simulating  quasi-steady  one-dimensional  laminar  premixed  flames  in  the 
configuration  of  a  jet  injected  in  air.  The  model  has  been  derived  in  the  framework  consistent  with  it  being 
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used  in  conjunction  with  a  reduced  chemical  kinetic  model  based  on  constituents  and  species.  As  such,  the 
formulation  includes  a  complete  mass-diffusion  matrix,  and  a  mixture  thermal  conductivity  computed  taking 
into  account  the  global  constituent  and  all  species.  Furthermore,  a  real-gas  equation  of  state  is  used. 

The  equations  are  solved  using  a  pseudo-unsteady  time  increment  procedure  with  initial  profiles  that  are 
a  guess  of  the  steady-state  profiles.  Results  have  been  presented  for  both  n- heptane/air  and  iso-octane/air 
premixed  flames  under  a  variety  of  initial  conditions.  Comparisons  with  experimental  results  showed  good 
agreement.  Further  simulations  are  necessary  to  explore  the  full  potential  of  the  model. 
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Appendix  A 

Miscellaneous  relationships  relevant  to  the  EOS  are 


dn 


EE  Clpn  (T) 

1  ^mix  —  ^  ^  ^p^p-} 


p  n  p 

where  indices  do  not  follow  here  the  Einstein  notation,  and 

dpn  —  (1  k  C^ppC^nm 

ctppiT)  =  0.457236(i?„Tc',p)^  1  -I-  Cp(l  —  \/Tred,p)  /Pc,p-, 

Cp  =  0.37464  -h  1.542260^  -  0.269920^, 

where  Tred,p  =  T/Tc^p,  Tc^p  is  the  critical  temperature  and  Op  is  the  acentric  factor.  Also, 

,RuTc,p 


bp  =  0.077796- 


PC,p 


(58) 

(59) 

(60) 
(61) 

(62) 


Tc,pn 


(1  kpn)  \/Tc^pRc,n  with  kpp  —  0, 

1/1/0  1  /o\  3 

^C,pn  “  g 

Tc  ,pn  Zc  ,pn 


-  _  (0,1/3  +  ,,1/3  V 

C,pn  —  2  iZc,p  T  Zc  n)  , 


PC.pn  — 


VCa 


(63) 

(64) 

(65) 

(66) 


with  Tred,pn  =  TjTc^pn,  Zc,p  being  the  critical  compression  factor  defined  as  Z  =  p/{pTRu/m),  vc,p  be¬ 
ing  the  critical  volume,  and  pc,p  being  the  critical  pressure,  kpn  is  an  empirical  mixing  parameter.  The 
relationship  between  the  parameters  kpn  and  fcp„  is 


(1  kpn)  —  (1  kpn) 


{vc,pvc,n)^^‘^ 

'^C,pn 


(67) 


Values  of  fcp„  are  listed  in  Table  3,  where  the  values  were  obtained  from  Reid  et  al.®  or  Knapp  et  al.i^ 
Otherwise,  for  pairs  involving  O2,  N2,  CO,  CO2  or  H2O,  ~  0. 
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Species 

m  (g/mol) 

Tc  (K) 

PC  (bar) 

vc  (cm® /mol) 

C7H16 

100.2 

540.2 

27.4 

432 

0.349 

CgHis 

114.23 

568.8 

24.9 

492 

0.396 

N2 

28.013 

126.26 

33.4 

89.8 

0.039 

H2O 

18.015 

647.3 

221 

57.1 

0.344 

CO2 

44.01 

304.1 

73.8 

93.9 

0.225 

O2 

32.0 

154.6 

50.43 

73.4 

0.025 

H 

1.008 

33.15 

(12) 

64.2 

(0.) 

H2 

2.016 

32.94 

12.84 

64.3 

-0.216 

CO 

28.01 

132.9 

35 

93.2 

0.066 

CH4 

16.04 

190.6 

45.94 

99.2 

0.0108 

H2O2 

34.015 

728 

220 

(77) 

(0.1) 

C2H2 

26.04 

308.3 

61.4 

112.7 

0.19 

C2H4 

28.054 

282.4 

50.4 

130.4 

0.089 

CH2O 

30.026 

408 

65.9 

(144) 

0.253 

Table  1.  Fuel  and  unsteady  light  species  properties. 


Radicals 

m  (g/mol) 

Tc  (K) 

PC  (bar) 

Vc  (cm® /mol) 

n 

0 

16.0 

116 

76 

35 

0. 

CH 

13.02 

183 

73 

58 

0. 

CH2 

14.027 

210 

66 

73 

0. 

CH3 

15.034 

221 

62 

82 

0. 

OH 

17.01 

167 

85 

45 

0. 

HCO 

29.02 

299 

70 

99 

0.25 

OOH 

33.01 

226 

83 

63 

0. 

HC2 

25.03 

291 

67 

100 

0.2 

C2H3 

27.044 

292 

57 

119 

0.1 

Table  2.  Quasi-steady  radicals  properties. 


11  of  14 


American  Institute  of  Aeronautics  and  Astronautics 


p 

n 

k' 

alkane 

alkane 

0.0 

alkane 

N2,  O2 

0.15 

alkane 

CO2 

0.11 

alkane 

H2O 

0.093-0.006nc 

alkane 

H2 

0.099nc-0.81 

H2 

N2,  O2 

0.12 

H2 

CO2 

-0.162 

CO2 

H2O 

0.095 

CO2 

N2,  O2 

-0.017 

H2O 

N2,  O2 

0.17 

Table  3. 


Values  of  k  for  species  pairs. 


ric  is  the  number  of  C  atoms  in  the  species. 
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200 


1900 


Figure  1.  (a)  Predictions  of  the  reduced  model  for  temperature  and  velocity,  (b)  corresponding  experimental 
data,  and  (c)  the  reactants’  mass  fractions.  Computations  and  experiments  are  for  heptane/air  with  0  =  1, 
po  =  1  bar.  To  =  350  K,  L  =  0.6  cm  and  (pu)o  =  0.18  g/cm^-s.  The  temperature  is  represented  by  the  solid  line, 
and  the  velocity  by  the  dashed  line. 
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Temperature, 


Figure  2.  Predictions  of  the  reduced  model  for  an  iso-octane/air  premixed  flame  under  the  following  conditions: 
^  p  =  10  bar,  To  =  300  K,  L  =  0.4  cm,  (pu)o  =  3.0  g/cm^-s.  (a)  Temperature  (solid  line)  and  velocity  (dashed 
line)  and  (b)  reactants  mass  fractions. 


Figure  3.  Predictions  of  the  reduced  model  for  an  heptane/air  premixed  flame  under  the  following  conditions: 
0  =  1,  p  =  10  bar.  To  =  500  K,  L  =  0.4  cm,  {pu)o  =  0.10  g/cm^-s.  (a)  Temperature  (solid  line)  and  velocity  (dashed 
line)  and  (b)  reactants  mass  fractions. 


Figure  4.  Predictions  of  the  reduced  model  for  an  heptane/air  premixed  flame  under  the  following  conditions: 
(p  =  1.9,  p  =  1  bar.  To  =  300  K,  L  =  0.5  cm,  (pu)o  =  0.0147  g/cm^-s.  (a)  Temperature  (solid  line)  and  velocity 
(dashed  line)  and  (b)  products  mass  fractions. 
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